#include<limits.h>
#include<math.h>
#include<iostream>
#include<fstream>
#include<sstream>
#include<cstdlib>
#include<vector>
#include<boost/foreach.hpp>
#include"mpi.h"

#ifdef USE_GPU
#include<CL/cl.hpp>
#include"../common/clsafe.h"
#endif
#include"hmm_impute_dimensions.h"
#include"../common/main.hpp"
#include"../common/analyzer.hpp"
#include"../common/io.hpp"
#include"../common/utility.hpp"
//#include"dimension.h"
#include "hmm_impute.hpp"

using namespace std;

void HMM_impute::debugfloat(string filename,float * arr,int len){
  ofstream ofs(filename.data());
  for(int i=0;i<len;++i){
    ofs<<i<<"\t"<<arr[i]<<endl;
    //ofs<<i<<":"<<arr[i]<<endl;
  }
  ofs.close();
}

void HMM_impute::flush2disk(ofstream & ofs_cache, packedshort_t * vec, int len){
  //cerr<<"Peeking before disk flush: ";
  //for(int i=len-10;i<len;++i) cerr<<" "<<vec[i].compressed[0]<<","<<vec[i].compressed[1];
  //cerr<<endl;
  char * charvec = reinterpret_cast<char *>(vec);
  //double start = clock();
  ofs_cache.write(charvec,sizeof(float)*len);
  //ofs<<"Disk flush: "<<(clock()-start)/CLOCKS_PER_SEC<<endl;
}

void HMM_impute::flush2ram(ifstream & ifs_cache, packedshort_t *  vec, int len,int row_offset){
  //double start = clock();
  unsigned long seekpos = row_offset*len*sizeof(packedshort_t);
  ofs<<"Seeking to "<<seekpos<<endl;
  ifs_cache.seekg(seekpos);
  //cerr<<"Reinterpreting cast\n";
  char * charvec = reinterpret_cast<char *>(vec);
  //cerr<<"Reading in data\n";
  ifs_cache.read(charvec,sizeof(packedshort_t)*len);
  //ofs<<"RAM flush: "<<(clock()-start)/CLOCKS_PER_SEC<<" at row "<<row_offset<<endl;
  vec = reinterpret_cast<packedshort_t *>(charvec);
  //cerr<<"Peeking after RAM flush: ";
  //for(int i=len-10;i<len;++i) cerr<<" "<<vec[i].compressed[0]<<","<<vec[i].compressed[1];
  //cerr<<endl;
}

void HMM_impute::flush2disk(ofstream & ofs_cache, float * vec, int len){
  //cerr<<"Peeking before disk flush: ";
  //for(int i=len-10;i<len;++i) cerr<<" "<<vec[i];
  //cerr<<endl;
  char * charvec = reinterpret_cast<char *>(vec);
  //double start = clock();
  ofs_cache.write(charvec,sizeof(float)*len);
  //ofs<<"Disk flush: "<<(clock()-start)/CLOCKS_PER_SEC<<endl;
}

void HMM_impute::flush2ram(ifstream & ifs_cache, float *  vec, int len,int row_offset){
  //double start = clock();
  unsigned long seekpos = row_offset*len*sizeof(float);
  ofs<<"Seeking to "<<seekpos<<endl;
  ifs_cache.seekg(seekpos);
  //cerr<<"Reinterpreting cast\n";
  char * charvec = reinterpret_cast<char *>(vec);
  //cerr<<"Reading in data\n";
  ifs_cache.read(charvec,sizeof(float)*len);
  //ofs<<"RAM flush: "<<(clock()-start)/CLOCKS_PER_SEC<<" at row "<<row_offset<<endl;
  vec = reinterpret_cast<float *>(charvec);
  //cerr<<"Peeking after RAM flush: ";
  //for(int i=len-10;i<len;++i) cerr<<" "<<vec[i];
  //cerr<<endl;
}

void HMM_impute::short2floatarr(packedshort_t * source,float * dest){
  //double  start = clock();
  int state =0;
  for(int cstate=0;cstate<compressed_states;++cstate){
    for(int i=0;i<2;++i){
      dest[state] = short2float(source[cstate].compressed[i]);
      ++state;
    }
  }
  //ofs<<"Converting short to float arr time: "<<(clock()-start)/CLOCKS_PER_SEC<<endl;
}

void HMM_impute::float2shortarr(float * source,packedshort_t * dest){
  //double  start = clock();
  int state =0;
  for(int cstate=0;cstate<compressed_states;++cstate){
    for(int i=0;i<2;++i){
      dest[cstate].compressed[i] = float2short(source[state]);
      ++state;
    }
  }
  //ofs<<"Converting float to short arr time: "<<(clock()-start)/CLOCKS_PER_SEC<<endl;
}

void HMM_impute::init(const ptree & pt){
  use_compressed = false;
  use_newfb = false;
  
  MPI_Init(0,NULL);
  MPI_Comm_size(MPI_COMM_WORLD,&mpi_numtasks);
  MPI_Comm_rank(MPI_COMM_WORLD,&mpi_rank);
  enabled_hosts = new bool[mpi_numtasks];
  if (mpi_rank){
    ostringstream oss;
    oss<<"debug.slave."<<mpi_rank;
    ofs.open(oss.str().data());
  }else{
    ostringstream oss;
    oss<<"debug.master";
    ofs.open(oss.str().data());
  }
  settings = new hmm_impute_settings_t;
  settings->marker_file = pt.get<string>("marker_file");
  settings->reads_file = pt.get<string>("reads_file");
  settings->cache_dir = pt.get<string>("cache_dir");
  settings->iterations = pt.get<int>("iterations");
  settings->max_block_size = pt.get<int>("max_block_size");
  settings->random_proportion = pt.get<float>("random_proportion");
  slaves = 0;
  int hostcounter = 0;
  BOOST_FOREACH(boost::property_tree::ptree::value_type host_object, pt.get_child("opencl_settings")){
    boost::property_tree::ptree host_tree = (boost::property_tree::ptree)host_object.second;
    int hostid = host_tree.get<int>("<xmlattr>.rank");
    bool enable_flag = !host_tree.get<string>("enabled").compare("true")?true:false;
    slaves+=(hostid>0 && enable_flag);
    enabled_hosts[hostcounter++] = enable_flag;
    if (hostid==mpi_rank){
      settings->host_enabled = enable_flag;
      ofs<<"This process has enable flag: "<<enable_flag<<endl;
      settings->use_gpu = !host_tree.get<string>("enable_opencl").compare("true")?true:false;
      settings->platform_id = host_tree.get<int>("platform_id");
      settings->device_id = host_tree.get<int>("device_id");
      ofs<<"This process wants device id: "<<settings->device_id<<endl;
      settings->kernel_path = host_tree.get<string>("kernel_path");
    }
  }
  ofs<<"There are a total of "<<slaves<<" slaves\n";
  if (settings->host_enabled){
    this->math = new MathUtils(123);
    ifstream ifs;
    ifs.open(settings->marker_file.data());
    if (!ifs.is_open()){
      cerr<<"Could not find marker map file: "<<settings->marker_file<<endl;
      exit(1);
    }
    string line;
    int lastpos=0;
    vector<string> rslist_temp;
    vector<int> distances_temp;
    while(getline(ifs,line)){
      istringstream iss(line);
      string chr,rs;
      int morgans,pos;
      iss>>chr>>rs>>morgans>>pos;
      rslist_temp.push_back(rs);
      int dist = lastpos==0?0:pos-lastpos;
      distances_temp.push_back(dist);
      lastpos = pos;
    }
    ifs.close();
    markers = rslist_temp.size();
    rslist = new string[markers];
    distances = new int[markers];
    for(int i=0;i<markers;++i){
      rslist[i] = rslist_temp[i];
      distances[i] = distances_temp[i];
    }
    
    ifs.open(settings->reads_file.data());
    if (!ifs.is_open()){
      cerr<<"Could not find seq read file: "<<settings->reads_file<<endl;
      exit(1);
    }
    vector<vector<float> > calls_temp;
    while(getline(ifs,line)){
      istringstream iss(line);
      vector<string> info(6);
      for(int i=0;i<6;++i) iss>>info[i];
      pedinfo.push_back(info);
      vector<float> callvec(markers*3);
      for(int i=0;i<markers*3;++i){
        iss>>callvec[i];
      }
      calls_temp.push_back(callvec);
    }
    ifs.close();
    //samples = 50;
    all_samples = calls_temp.size();
    sample_mask = new bool[all_samples];
    random_samples = static_cast<int>(settings->random_proportion*all_samples) ;
    ofs<<"Random samples for template haplotypes: "<<random_samples<<endl;
    seq_error = .002;
    base_error_rate = .01;
    base_xover_rate = .01;
    site_errors = new float[markers];
    site_thetas = new float[markers];
    for(int j=0;j<markers;++j) {
      site_errors[j] = base_error_rate;
      site_thetas[j] = base_xover_rate;
    }
    // to cache the computation of binomial distribution probailities
    // 3 possible true genotypes, 16 possible counts for the A and B allele
    // for G=2,A=3,B=5 we dereference as 2*256+3*16+5;
    binomial_matrix = new float[3*16*16];
    penetrance_matrix = new float[9*markers];
    for(int a=0;a<16;++a){
      for(int b=0;b<16;++b){
        int c = a+b;
        binomial_matrix[0*256+a*16+b] = binom_prob(c,a,seq_error);
        binomial_matrix[1*256+a*16+b] = binom_prob(c,a,.5);
        binomial_matrix[2*256+a*16+b] = binom_prob(c,b,seq_error);
      }
    }
    call_probs = new float[all_samples * markers*3]; // subject major
    dose_matrix = new int[all_samples * markers];
    for(int i=0;i<all_samples;++i){
      for(int j=0;j<markers;++j){
        // there are three call probs per marker
        call_probs[i*markers*3+j*3] = calls_temp[i][j*3];
        call_probs[i*markers*3+j*3+1] = calls_temp[i][j*3+1];
        call_probs[i*markers*3+j*3+2] = calls_temp[i][j*3+2];
        dose_matrix[i*markers+j] = 0;
      }
    }
    ofs<<"There are a total of "<<markers<<" markers across "<<all_samples<<" subjects.\n";
    ofs<<"Input files loaded.\n";
    total_haplotypes = 2*all_samples;
    total_template_haplotypes = 2*(random_samples-1);
    total_template_haplotype_pairs = total_template_haplotypes/2;
    states = (total_template_haplotypes)*(total_template_haplotypes);  // cannot copy from self
    compressed_states = states/2;
    int state_vector_bytes = states * 4;
    matrix_compression_ratio = settings->max_block_size/state_vector_bytes;
    ofs<<"To support block sizes of "<<(settings->max_block_size/1e6)<<" MB, matrix will have "<<matrix_compression_ratio<<" rows.\n";
    state_offsets = new int[total_template_haplotypes];
    for(int i=0;i<total_template_haplotypes;++i){
      state_offsets[i] = i*total_template_haplotypes;
    }
    state_haplo_indices = new int[states*2];
    ofs<<"There are "<<states<<" states.\n";
    ofs<<"Allocated state haplo indices of size: "<<states*2<<endl;
    int state = 0;
    for(int i=0;i<total_template_haplotypes;++i){
      for(int j=0;j<total_template_haplotypes;++j){
        state_haplo_indices[state*2] = i;
        state_haplo_indices[state*2+1] = j;
        ++state;
      }
    }
    // for odd vs even markers
    haplotype_marginals = new float[total_template_haplotypes * 2];
    estimated_haplotype_matrix = new int[markers * total_haplotypes]; 
    estimated_haplotype_matrix_subjectmajor = new int[total_haplotypes*markers]; 
    short_template_haplotypes = new packedshort_t[markers * total_template_haplotype_pairs]; 
    //template_haplotypes = new int[markers * total_template_haplotypes]; 
    init_haplotypes();
    //estimated_haplotype_pair = new int[markers*2]; 
    //for(int i=0;i<2*markers;++i) estimated_haplotype_pair[i] = 0;
    //setup_haplotypes();
    //trans_mat = new float[states* states];
    site_matches = new int[markers];
    site_mismatches = new int[markers];
    site_xovers = new int[markers];
    geno_probs = new float[all_samples * 3 * markers];
    
    current_block_id = -1;
    likelihood_log = -9999;
    //forward_rescaling = new int[markers];
    //backward_rescaling = new int[markers];
    matrix_signposts = markers/matrix_compression_ratio + (markers%matrix_compression_ratio!=0); 
      //emission_mat_compressed = new float[matrix_signposts*states];
    temp_prob_vec = new float[states];
    emission_vec = new float[states];
    state = 0;
    if(use_compressed){
      short_previous_prob_vec = new packedshort_t[compressed_states];
      short_current_prob_vec = new packedshort_t[compressed_states];
      short_prior_vec = new packedshort_t[compressed_states];
      for(int cstate = 0;cstate<compressed_states;++cstate){
        for(int i=0;i<2;++i){
          short_prior_vec[cstate].compressed[i] = state_haplo_indices[state*2]==state_haplo_indices[state*2+1]? float2short(1./(total_haplotypes*total_haplotypes)):float2short(1);
          ++state;
        }
      }
      short_temp_prob_vec = new packedshort_t[states];
      short_backward_mat_block = new packedshort_t[matrix_compression_ratio*compressed_states];
    }else{
      previous_prob_vec = new float[states];
      current_prob_vec = new float[states];
      prior_vec = new float[states];
      for(int state=0;state<states;++state){
        prior_vec[state] = state_haplo_indices[state*2]==state_haplo_indices[state*2+1]? 1./(total_haplotypes*total_haplotypes):1.;
      }
      backward_mat_block = new float[matrix_compression_ratio*states];
    }
    
    //emission_mat_block = new float [matrix_compression_ratio*states];
    if(settings->use_gpu){
      init_gpu();
    }
    int rc;
    MPI_Type_contiguous(all_samples*markers,MPI_INT,&MatrixIntArrayType);
    MPI_Type_commit(&MatrixIntArrayType);
    MPI_Type_contiguous(markers,MPI_INT,&VectorIntArrayType);
    MPI_Type_commit(&VectorIntArrayType);
    MPI_Type_contiguous(markers,MPI_FLOAT,&VectorFloatArrayType);
    MPI_Type_commit(&VectorFloatArrayType);
    if (mpi_rank==0){
      samples_per_slave = new int[slaves];
      int division = all_samples/slaves;
      int offset = 0;
      int slave=0;
      for(int dest=1;dest<mpi_numtasks;++dest){
        if (enabled_hosts[dest]){
          int tasks = (slave<slaves-1)?division:(all_samples-(slaves-1)*division);
          samples_per_slave[slave] = tasks;
          cerr<<"Slave "<<slave<<" has "<<tasks<<" tasks.\n";
          rc = MPI_Send(&offset,1,MPI_INT,dest,TAG_INIT_SUBJECT_OFFSET,MPI_COMM_WORLD);
          rc = MPI_Send(&tasks,1,MPI_INT,dest,TAG_INIT_SAMPLESIZE,MPI_COMM_WORLD);
          offset += tasks;
          ++slave;
        }
      }
    }else{
      for(int i=0;i<matrix_compression_ratio;++i){
        for(int j=0;j<states;++j){
          backward_mat_block[i*states+j] = 0;
        }
        for(int j=0;j<compressed_states;++j){
          for(int k=0;k<2;++k){
            //short_backward_mat_block[i*compressed_states+j].compressed[k] = 0;
          }
        }
      }
      //for(int i=0;i<matrix_compression_ratio*compressed_states;++i) backward_mat_block[i] = 0;
      rc = MPI_Recv(&slave_sample_offset,1,MPI_INT,0,TAG_INIT_SUBJECT_OFFSET,MPI_COMM_WORLD,&stat);
      rc = MPI_Recv(&slave_sample_size,1,MPI_INT,0,TAG_INIT_SAMPLESIZE,MPI_COMM_WORLD,&stat);
      cerr<<"I am slave "<<mpi_rank<<" and have offset: "<<slave_sample_offset<<" and samplesize: "<<slave_sample_size<<endl;
      ostringstream oss;
      oss<<settings->cache_dir<<"/rank."<<mpi_rank<<".bin";
      cachefilename = oss.str();
    } // end if master or slave
  }
  ofs<<"Done with init()\n";
}

HMM_impute::HMM_impute(){
}

void HMM_impute::init_haplotypes(){
  for(int j=0;j<markers;++j){
    int haploindex = 0;
    for(int i=0;i<all_samples;++i){
      float r = math->RandomUniform();
      float p1 = call_probs[i*3*markers+3*j];
      float p2 = call_probs[i*3*markers+3*j+1];
      //float p3 = call_probs[i*3*markers+3*j+2];
      if (r<p1){
        estimated_haplotype_matrix[j*total_haplotypes+haploindex] = 0;
        estimated_haplotype_matrix[j*total_haplotypes+haploindex+1] = 0;
      }else if (r<p1+p2){
        bool bit = math->RandomUniform()>.5;
        estimated_haplotype_matrix[j*total_haplotypes+haploindex] = bit;
        estimated_haplotype_matrix[j*total_haplotypes+haploindex+1] = !bit;
      }else{
        estimated_haplotype_matrix[j*total_haplotypes+haploindex] = 1;
        estimated_haplotype_matrix[j*total_haplotypes+haploindex+1] = 1;
      }
      haploindex+=2;
    }
  }
  ofs<<"Transposing haplotypes to MPI communicable data structure\n";
  for(int i=0;i<total_haplotypes;++i){
    for(int j=0;j<markers;++j){
      estimated_haplotype_matrix_subjectmajor[i*markers+j] = estimated_haplotype_matrix[j*total_haplotypes+i];
    }
  }
  bool debug_haplo = false;
  if (debug_haplo){
    ofstream ofs("debug.haplo");
    for(int j=0;j<markers;++j){
      for(int i=0;i<total_haplotypes;++i){
        cerr<< estimated_haplotype_matrix[j*total_haplotypes+i];
      }
      cerr<<endl;
    }
    ofs.close();
  }
}

void HMM_impute::init_iteration(){
  // normally the MPI master will broadcast the best pool of estimated haplotypes to the original matrix.
  int rc;
  if (mpi_rank==0){
    cerr<<"Transferring best haplotype estimates and site thetas to the slaves\n";
    for(int dest=1;dest<mpi_numtasks;++dest){
      if (enabled_hosts[dest]){
        rc = MPI_Send(estimated_haplotype_matrix_subjectmajor,1,MatrixIntArrayType,dest,TAG_NEW_HAPLOTYPE_POOL,MPI_COMM_WORLD);
        rc = MPI_Send(site_thetas,1,VectorFloatArrayType,dest,TAG_NEW_SITE_RATES,MPI_COMM_WORLD);
        rc = MPI_Send(site_errors,1,VectorFloatArrayType,dest,TAG_NEW_SITE_RATES,MPI_COMM_WORLD);
      }
    }
  }else{
    ofs<<"Receiving best haplotype and theta estimates to the haplotype source pool\n";
    rc = MPI_Recv(estimated_haplotype_matrix_subjectmajor,all_samples*markers,MPI_INT,0,TAG_NEW_HAPLOTYPE_POOL,MPI_COMM_WORLD,&stat);
    ofs<<"Transposing haplotypes from MPI communicable data structure\n";
    for(int i=0;i<total_haplotypes;++i){
      for(int j=0;j<markers;++j){
        estimated_haplotype_matrix[j*total_haplotypes+i] = estimated_haplotype_matrix_subjectmajor[i*markers+j];
     }
    }
//    for(int i=0;i<samples;++i){
//      ofs<<"subject "<<i<<":";
//      for(int j=0;j<10;++j){
//        ofs<<estimated_haplotype_matrix_subjectmajor[(i*2)*markers+j];
//      }
//      ofs<<endl;
//    }
//    ofs.close();
//    exit(0);
    rc = MPI_Recv(site_thetas,markers,MPI_FLOAT,0,TAG_NEW_SITE_RATES,MPI_COMM_WORLD,&stat);
    rc = MPI_Recv(site_errors,markers,MPI_FLOAT,0,TAG_NEW_SITE_RATES,MPI_COMM_WORLD,&stat);
  // update the penetrance matrix based on the last site of error and theta estimates
    for(int marker=0;marker<markers;++marker){
      // both homozygote and matching
      penetrance_matrix[9*marker+0*3+0] = penetrance_matrix[9*marker+2*3+2] = (1-site_errors[marker])*(1-site_errors[marker]);
      // G homozygote but S is het
      penetrance_matrix[9*marker+0*3+1] = penetrance_matrix[9*marker+2*3+1] = (1-site_errors[marker])*(site_errors[marker]);
      // both G homozygotes and opposite to S
      penetrance_matrix[9*marker+0*3+2] = penetrance_matrix[9*marker+2*3+0] = (site_errors[marker])*(site_errors[marker]);
      // G het and both S homos
      penetrance_matrix[9*marker+1*3+0] = penetrance_matrix[9*marker+1*3+2] = 2*(site_errors[marker])*(1-site_errors[marker]);
      // both het and matching
      penetrance_matrix[9*marker+1*3+1] = (1-site_errors[marker])*(1-site_errors[marker])+(site_errors[marker])*(site_errors[marker]);
    }
    if (settings->use_gpu){
#ifdef USE_GPU
      err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.penetrance_matrix_buffer, CL_TRUE, 0,  sizeof(float)*(markers*9), penetrance_matrix , NULL, NULL );
      clSafe(err, "write buffer for penetrance matrix\n");
      err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.site_thetas_buffer, CL_TRUE, 0,  sizeof(float)*(markers), site_thetas, NULL, NULL );
      clSafe(err, "write buffer for site_thetas\n");
#endif
    }
  }
  // reset the number of error and cross over tallies for this round
  for(int j=0;j<markers;++j){
    site_matches[j] = 0;
    site_mismatches[j] = 0;
    site_xovers[j] = 0;
  }
}

void HMM_impute::finalize_subject(){
  ofs<<"Updating doses for subject "<<current_subject<<endl;
  for(int marker=0;marker<markers;++marker){
    dose_matrix[current_subject*markers+marker] += estimated_haplotype_matrix[marker*total_haplotypes+current_subject*2] + estimated_haplotype_matrix[marker*total_haplotypes+current_subject*2+1];
    if (marker>markers-50){
      ofs<<dose_matrix[current_subject*markers+marker];
    }
  }
  ofs<<endl;
}

void HMM_impute::write_dosage(){
  ostringstream oss;
  oss<<"results.dosages."<<iter;
  ofstream ofs_dose(oss.str().data());
  char dosestr[5];
  for(int i=0;i<all_samples;++i){
    for(int j=0;j<6;++j){
      if (j) ofs_dose<<" ";
      ofs_dose<<pedinfo[i][j];
    }
    for(int marker=0;marker<markers;++marker){
      //ofs_dose<<geno_probs[i*3*markers+3*marker+1]+2*geno_probs[i*3*markers+3*marker+2];
      //ofs_dose<< dose_matrix[i*markers+marker];
     //ofs_dose<<(2*(iter+1)-dose_matrix[i*markers+marker])/(iter+1);
     float d = abs(1.*(2*(iter+1)-dose_matrix[i*markers+marker])/(iter+1)-2.);
     sprintf(dosestr,"%1.3f",d);
     ofs_dose<<" "<<dosestr;
     //ofs_dose<<" "<<dose_matrix[i*markers+marker]<<","<<dosestr;
     //ofs_dose<<dose_matrix[i*markers+marker];
    }
    ofs_dose<<endl;
  }
  ofs_dose.close(); 
}

void HMM_impute::finalize(){
}

void HMM_impute::finalize_iteration(){
  int rc,rpc_code;
  if (mpi_rank==0){
    rpc_code = RPC_REQUEST_HAPLO;
    int offset = 0;
    int slave = 0; 
    for(int dest=1;dest<mpi_numtasks;++dest){
      if (enabled_hosts[dest]){
        rc = MPI_Send(&rpc_code,1,MPI_INT,dest,RPC,MPI_COMM_WORLD);
        for(int i = offset;i<offset+samples_per_slave[slave];++i){
          ofs<<"Receiving updated haplotypes for subject "<<i<<endl;
          rc = MPI_Recv(estimated_haplotype_matrix_subjectmajor+(2*i*markers),markers,MPI_INT,dest,TAG_NEW_HAPLOTYPE,MPI_COMM_WORLD,&stat);
          rc = MPI_Recv(estimated_haplotype_matrix_subjectmajor+((2*i+1)*markers),markers,MPI_INT,dest,TAG_NEW_HAPLOTYPE,MPI_COMM_WORLD,&stat);
          for(int j=0;j<50;++j){
           ofs<<estimated_haplotype_matrix_subjectmajor[2*i*markers+j]+estimated_haplotype_matrix_subjectmajor[(2*i+1)*markers+j];
          }
          ofs<<endl;
        }
        offset+=samples_per_slave[slave];
        ++slave;
      }
    }
    // update the site specific error and crossover rates
    //write_geno_posterior();
  }else{
    ofs<<"Transposing haplotypes to MPI communicable data structure\n";
    for(int i=0;i<total_haplotypes;++i){
      for(int j=0;j<markers;++j){
        estimated_haplotype_matrix_subjectmajor[i*markers+j] = estimated_haplotype_matrix[j*total_haplotypes+i];
      }
    }
    rc = MPI_Recv(&rpc_code,1,MPI_INT,0,RPC,MPI_COMM_WORLD,&stat);
    if (rpc_code==RPC_REQUEST_HAPLO){
      ofs<<"Sending updated haplotypes\n";
      for(int i=slave_sample_offset;i<slave_sample_offset+slave_sample_size;
      ++i){
        ofs<<"Sending updated haplotypes for subject "<<i<<endl;
        rc = MPI_Send(estimated_haplotype_matrix_subjectmajor+(2*i*markers),1,VectorIntArrayType,0,TAG_NEW_HAPLOTYPE,MPI_COMM_WORLD);
        rc = MPI_Send(estimated_haplotype_matrix_subjectmajor+((2*i+1)*markers),1,VectorIntArrayType,0,TAG_NEW_HAPLOTYPE,MPI_COMM_WORLD);
        for(int j=0;j<50;++j){
         ofs<<estimated_haplotype_matrix_subjectmajor[2*i*markers+j]+estimated_haplotype_matrix_subjectmajor[(2*i+1)*markers+j];
        }
        ofs<<endl;
      }
    }else{
      cerr<<"Something is wrong with the request code, expecting haplotype update\n";
      exit(0);
    }
  }
  //int rc,rpc_code;
  if (mpi_rank==0){
    rpc_code = RPC_REQUEST_DOSE;
    int offset = 0;
    int slave = 0;
    for(int dest = 1;dest<mpi_numtasks;++dest){
      if (enabled_hosts[dest]){
        rc = MPI_Send(&rpc_code,1,MPI_INT,dest,RPC,MPI_COMM_WORLD);
        for(int i = offset;i<offset+samples_per_slave[slave];++i){
          ofs<<"Receiving updated doses for subject "<<i<<endl;
          rc = MPI_Recv(dose_matrix+(i*markers),markers,MPI_INT,dest,TAG_NEW_DOSE,MPI_COMM_WORLD,&stat);
          for(int j=markers-50;j<markers;++j){
            ofs<<dose_matrix[i*markers+j];
          }
          ofs<<endl;
        }
        offset+=samples_per_slave[slave];
        ++slave; 
      }
    }
    write_dosage();
  }else{
    rc = MPI_Recv(&rpc_code,1,MPI_INT,0,RPC,MPI_COMM_WORLD,&stat);
    if (rpc_code==RPC_REQUEST_DOSE){
      for(int i=slave_sample_offset;i<slave_sample_offset+slave_sample_size;
      ++i){
        ofs<<"Sending updated doses for subject "<<i<<endl;
        rc = MPI_Send(dose_matrix+(i*markers),1,VectorIntArrayType,0,TAG_NEW_DOSE,MPI_COMM_WORLD);
        for(int j=markers-50;j<markers;++j){
          ofs<<dose_matrix[i*markers+j];
        }
        ofs<<endl;
      }
    }else{
      cerr<<"Something is wrong with the request code, expecting dose update\n";
      exit(0);
    }
  } 
  if (mpi_rank==0){
    rpc_code = RPC_REQUEST_SITE_RATES;
    int temp_xovers[markers];
    int temp_matches[markers];
    int temp_mismatches[markers];
    for(int dest = 1;dest<mpi_numtasks;++dest){
      if (enabled_hosts[dest]){
        rc = MPI_Send(&rpc_code,1,MPI_INT,dest,RPC,MPI_COMM_WORLD);
        ofs<<"Receiving updated xovers from rank "<<dest<<endl;
        rc = MPI_Recv(temp_xovers,markers,MPI_INT,dest,TAG_NEW_SITE_RATES,
        MPI_COMM_WORLD,&stat);
        ofs<<"Receiving updated matches from rank "<<dest<<endl;
        rc = MPI_Recv(temp_matches,markers,MPI_INT,dest,TAG_NEW_SITE_RATES,
        MPI_COMM_WORLD,&stat);
        ofs<<"Receiving updated mismatches from rank "<<dest<<endl;
        rc = MPI_Recv(temp_mismatches,markers,MPI_INT,dest,TAG_NEW_SITE_RATES,
        MPI_COMM_WORLD,&stat);
        for(int j=0;j<markers;++j) {
          site_xovers[j]+=temp_xovers[j];
          site_matches[j]+=temp_matches[j];
          site_mismatches[j]+=temp_mismatches[j];
          //ofs<<" "<<j<<","<<site_xovers[j];
        }
        //ofs<<endl;
      }
    }
    update_rates();
  }else{
    rc = MPI_Recv(&rpc_code,1,MPI_INT,0,RPC,MPI_COMM_WORLD,&stat);
    if (rpc_code==RPC_REQUEST_SITE_RATES){
      ofs<<"Sending updated thetas,matches,mismatches"<<endl;
      rc = MPI_Send(site_xovers,1,VectorIntArrayType,0,TAG_NEW_SITE_RATES,MPI_COMM_WORLD);
      rc = MPI_Send(site_matches,1,VectorIntArrayType,0,TAG_NEW_SITE_RATES,MPI_COMM_WORLD);
      rc = MPI_Send(site_mismatches,1,VectorIntArrayType,0,TAG_NEW_SITE_RATES,MPI_COMM_WORLD);
    }else{
      cerr<<"Something is wrong with the request, expecting theta update\n";
      exit(0);
    }
  } 
  // end if mpi master or slave
}


void HMM_impute::write_geno_posterior(){
  ostringstream oss;
  oss<<"results.posteriors."<<iter;
  ofstream ofs(oss.str().data());
  for(int i=0;i<all_samples;++i){
    for(int j=0;j<6;++j){
      if (j) cerr<<" ";
      cerr<<pedinfo[i][j];
    }
    for(int marker=0;marker<markers;++marker){
      cerr<<" ";
      //int g = estimated_haplotype_matrix[(i*2)*markers+j] +  estimated_haplotype_matrix[(i*2+1)*markers+j];
      //cerr<<"\t"<<g;
      cerr<<geno_probs[i*3*markers+3*marker] <<","<<geno_probs[i*3*markers+3*marker+1] <<","<<geno_probs[i*3*markers+3*marker+2];
    }
    cerr<<endl;
  }
  ofs.close();
}

void HMM_impute::init_subject(){
  bool random = settings->random_proportion<1;
  if (random){
    int random_count = 0;
    for(int i=0;i<all_samples;++i) sample_mask[i] = false;
    ofs<<"Sampling random subjects...\n";
    while(random_count<(random_samples-1)){
      int randindex = static_cast<int>(math->RandomUniform()*all_samples);
      if (randindex<all_samples && randindex!=current_subject && 
      !sample_mask[randindex]){
        sample_mask[randindex] = true;
        ++random_count;
      }
    }
    ofs<<"Done sampling random subjects.\n";
  }else{
    ofs<<"In debug mode: using all subjects\n";
    for(int i=0;i<all_samples;++i){
      sample_mask[i] = i!=current_subject;
    }
  }
  trace = current_subject==-20;
  for(int marker=0;marker<markers;++marker){
    int happairindex = 0;
    for(int subject = 0;subject<all_samples;++subject){
      if (sample_mask[subject]){
        //cerr<<"Bsubject:"<<subject<<":"<<estimated_haplotype_matrix[marker*total_haplotypes+(subject*2)]<<original_haplotype_matrix[marker*total_haplotypes+(subject*2)]<<estimated_haplotype_matrix[marker*total_haplotypes+(subject*2+1)]<<original_haplotype_matrix[marker*total_haplotypes+(subject*2+1)]<<endl;
        short_template_haplotypes[marker*total_template_haplotype_pairs+happairindex].compressed[0] = estimated_haplotype_matrix[marker*total_haplotypes+(subject*2)];
        short_template_haplotypes[marker*total_template_haplotype_pairs+happairindex].compressed[1] = estimated_haplotype_matrix[marker*total_haplotypes+(subject*2+1)];
        //template_haplotypes[marker*total_template_haplotypes+hapindex] = original_haplotype_matrix[marker*total_haplotypes+(subject*2)];
        //template_haplotypes[marker*total_template_haplotypes+(hapindex+1)] = original_haplotype_matrix[marker*total_haplotypes+(subject*2+1)];
        ++happairindex;
      }
    }
  }
  //ofstream ofs("debug.haplo.2");
  //cerr<<"Current Subject: "<<current_subject<<endl;
  //for(int j=0;j<markers;++j){
    //int hapindex = 0;
    //for(int i=0;i<samples;++i){
      //if (current_subject==i){
        //cerr<<"??";
      //}else{
        //cerr<<template_haplotypes[j*total_haplotypes+hapindex]<<template_haplotypes[j*total_haplotypes+hapindex+1];
        //hapindex+=2;
      //}
    //}
    //cerr<<endl;
  //}
  //ofs.close();
  //exit(0);
  // make sure we load a new  posterior block
  current_block_id = -1;
  if (settings->use_gpu){
  #ifdef USE_GPU
    err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.current_subject_buffer, CL_TRUE, 0,  sizeof(int)*(1), &current_subject , NULL, NULL );
    clSafe(err, "write buffer for current_person\n");
    err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.short_template_haplotypes_buffer, CL_TRUE, 0,  sizeof(packedshort_t)*(total_template_haplotype_pairs*markers), short_template_haplotypes , NULL, NULL );
    clSafe(err, "write buffer for haplotypes\n");
  #endif
  }
}
  
void HMM_impute::run(){
  if (settings->host_enabled){
    for(iter = 0;iter<settings->iterations;++iter){
      init_iteration();
      if (mpi_rank>0){
        for(current_subject = slave_sample_offset;current_subject<slave_sample_offset+slave_sample_size;++current_subject){
          ofs<<"Running subject: "<<current_subject<<endl;
          double start1 = clock();
          init_subject();
          double start = clock();
          compute_matrix_signposts();
          ofs<<"Time for signposts: "<<(clock()-start)/CLOCKS_PER_SEC<<endl;
          start = clock();
          sample_chromosomes();
          ofs<<"Time for sample chrs: "<<(clock()-start)/CLOCKS_PER_SEC<<endl;
          finalize_subject();
          double end = (clock()-start1)/CLOCKS_PER_SEC;
          ofs<<"Time for subject "<<current_subject<<": "<<end<<endl;
          cerr<<"Time for subject "<<current_subject<<": "<<end<<endl;
        } // loop across subjects
      }else{
        cerr<<"Running iteration: "<<iter<<endl;
      }
      finalize_iteration();
    }// loop iterations
    --iter;
    finalize();
  }
  ofs<<"Finalizing MPI\n";
  MPI_Finalize();
  ofs.close();
}


HMM_impute::~HMM_impute(){
  cerr<<"Destructor\n";
}

float HMM_impute::binom_prob(int n,int k,float p){
  float num = 1;
  for(int i=n;i>(n-k);--i){
    num*=i;
  }
  float den=1;
  for(int i=k;i>1;--i){
    den*=i;
  }
  float binom = num/den*pow(p,k)*pow(1.-p,n-k);
  return binom;
}

void HMM_impute::update_rates(){
    unsigned int base_length = 0;
    unsigned int base_xovers = 1;
    unsigned int base_mismatches = 0;
    unsigned int base_matches = 0;
    ostringstream oss_theta;
    ostringstream oss_error;
    oss_theta<<"theta."<<iter;
    oss_error<<"error."<<iter;
    ofstream ofs_theta(oss_theta.str().data());
    ofstream ofs_error(oss_error.str().data());

    float scale = 1./total_haplotypes;
    for(int j=0;j<markers;++j){
      //cerr<<"marker "<<j<<" match/mismatch: "<<site_matches[j]<<"/"<<site_mismatches[j]<<endl;
      if (site_mismatches[j]>2){
        site_errors[j] = 1.*site_mismatches[j]/(site_mismatches[j]+site_matches[j]);
        //if (site_errors[j]<min_error_rate) site_errors[j] = min_error_rate;
        ofs_error<<j<<"\t"<<site_errors[j]<<endl;
        //cerr<<j<<"\t"<<site_errors[j]<<endl;
        //cerrcerr<<"Setting site error at "<<j<<" to "<<site_errors[j]<<endl;
      }else{
        base_matches += site_matches[j];
        base_mismatches += site_mismatches[j];
      }
      if (site_xovers[j]<=1){
        base_xovers+=site_xovers[j];
        ++base_length;
        //base_length+=distances[j];
      }
    }
    // update the global error and crossover rates
    if (base_length<1) base_length = 1;
    //cerr<<"Base xovers and len: "<<base_xovers<<","<<base_length<<endl;
    cerr<<"Scale is "<<scale<<endl;
    base_xover_rate = base_xovers * scale / base_length;
    //if (base_xover_rate<min_xover_rate) base_xover_rate=min_xover_rate;
    cerr<<"base_length: "<<base_length<<" Base xover rate is now "<<base_xover_rate<<endl;
    int d = base_mismatches+base_matches;
    //if (d<1) d=1;
    base_error_rate = 1.*base_mismatches/d ;
    //if (base_error_rate<min_error_rate) base_error_rate=min_error_rate;
    cerr<<"Base error rate is now "<<base_error_rate<<endl;
    for(int j=0;j<markers;++j){
      if (site_mismatches[j]<=2){
        site_errors[j] = base_error_rate; 
        ofs_error<<j<<"\t"<<site_errors[j]<<endl;
      }
      if (site_xovers[j]<=1){
        //float haldane = .5*(1-exp(-2 * distances[j]*1e-8));
        site_thetas[j] =  base_xover_rate;
        //site_thetas[j] =  haldane;
      }else{
        site_thetas[j] =  site_xovers[j]*scale;
      }
      ofs_theta<<j<<"\t"<<site_thetas[j]<<endl;
        //cerr<<"Setting "<<j<<" : "<<site_thetas[j]<<" with Haldane: "<<haldane<<endl;
    } 
    ofs_theta.close();
    ofs_error.close();
}

void HMM_impute::compute_emission_vec(int marker, float * emission_vec){
//  float * calls = call_probs + (current_subject*markers*3+marker*3);
//  float * pens = penetrance_matrix + (marker*9);
//  int * haps = template_haplotypes + (marker*total_template_haplotypes);
  float calls[3];
  for(int i=0;i<3;++i) {
    calls[i] = call_probs[current_subject*markers*3+marker*3+i];
  }
  float pens[9];
  for(int i=0;i<9;++i) {
    pens[i] = penetrance_matrix[marker*9+i];
  }
  short int haps[total_template_haplotypes];
  for(int i=0;i<total_template_haplotype_pairs;++i) {
    for(int j=0;j<2;++j) {
      haps[i*2+j] = short_template_haplotypes[marker*total_template_haplotype_pairs+i].compressed[j];
    }
  }
  for(int hap1=0;hap1<total_template_haplotypes;++hap1){
    int state_offset = hap1*total_template_haplotypes;
    int s = haps[hap1]+haps[hap1];
    //emission_vec[state_offset+hap1]=s;
    emission_vec[state_offset+hap1]=calls[0] * pens[s] + calls[1] * pens[3+s] + calls[2] * pens[6+s];
    for(int hap2=0;hap2<hap1;++hap2){
      int s = haps[hap1]+haps[hap2];
      //emission_vec[state_offset+hap2]=s;
      emission_vec[state_offset+hap2]=calls[0] * pens[s] + calls[1] * pens[3+s] + calls[2] * pens[6+s];
      int otherstate = hap2*total_template_haplotypes+hap1;
      //emission_vec[otherstate] = 999;
      emission_vec[otherstate] = emission_vec[state_offset+hap2];

    }
  }
}


void HMM_impute::impute_and_evaluate(int marker,int hap1, int hap2){
  //cerr<<"Subject "<<current_subject<<" At marker: "<<marker<<endl;
  //cerr<<"haps "<<hap1<<","<<hap2<<endl;
  //int h1 = template_haplotypes[marker*total_template_haplotypes+hap1];
  //int h2 = template_haplotypes[marker*total_template_haplotypes+hap2];
  short int h1 = short_template_haplotypes[marker*total_template_haplotype_pairs+(hap1/2)].compressed[hap1%2];
  short int h2 = short_template_haplotypes[marker*total_template_haplotype_pairs+(hap2/2)].compressed[hap2%2];
  //cerr<<"h1/h2 "<<h1<<"/"<<h2<<endl;
  int s = h1+h2;
  //cerr<<"S is "<<s<<endl;
  //float pen0 = call_probs[current_subject*markers*3+marker*3];
  float pen0 = penetrance_matrix[9*marker+0*3+s] * call_probs[current_subject*markers*3+marker*3];
  //float pen1 = call_probs[current_subject*markers*3+marker*3+1];
  float pen1 = penetrance_matrix[9*marker+1*3+s] * call_probs[current_subject*markers*3+marker*3+1];
  //float pen2 =  call_probs[current_subject*markers*3+marker*3+2];
  float pen2 = penetrance_matrix[9*marker+2*3+s] * call_probs[current_subject*markers*3+marker*3+2];
  //cerr<<"Pens: "<<pen0<<","<<pen1<<","<<pen2<<endl;
  float normalizer = pen0+pen1+pen2;
  float r = math->RandomUniform() * normalizer;
  if (r<pen0){
      //cerr<<"Imputing 11\n";
    //cerr<<"chose a homo wt\n";
    estimated_haplotype_matrix[marker*total_haplotypes+current_subject*2] = 0;
    estimated_haplotype_matrix[marker*total_haplotypes+current_subject*2+1] = 0;
  }else if (r<pen0+pen2){
      //cerr<<"Imputing 22\n";
    //cerr<<"chose a homo mutant\n";
    estimated_haplotype_matrix[marker*total_haplotypes+current_subject*2] = 1;
    estimated_haplotype_matrix[marker*total_haplotypes+current_subject*2+1] = 1;
  }else{ // the estimated haplotype should be a het  
    //cerr<<"chose a het\n";
    //cerr<<"het for "<<current_subject<<" at marker "<<marker<<endl;
    float r = math->RandomUniform();
    if ( h1 != h2 ){ // if template is also a het
      float error = pow(site_errors[marker],2)/(pow(site_errors[marker],2)+pow(1.-site_errors[marker],2)) ;
      //if (error>r) cerr<<"Swapping het\n";
      //else cerr<<"Retaining het\n";
      estimated_haplotype_matrix[marker*total_haplotypes+current_subject*2] = error>r?!h1:h1;
      estimated_haplotype_matrix[marker*total_haplotypes+current_subject*2+1] = error>r?!h2:h2;
    }else{ // the template is homo, so make a random het
      //cerr<<"Imputing het\n";
      bool b1 = r>.5;
      bool b2 = !b1;
      estimated_haplotype_matrix[marker*total_haplotypes+current_subject*2] = b1;
      estimated_haplotype_matrix[marker*total_haplotypes+current_subject*2+1] = b2;
    }
    //cerr<<"Made "<<estimated_haplotype_pair[2*marker]<<" and "<<estimated_haplotype_pair[2*marker+1]<<endl;
  }
  //cerr<<"marker "<<marker<<endl;
  // compare the differences between the imputed and the templates
  int diff = abs(estimated_haplotype_matrix[marker*total_haplotypes+current_subject*2]-h1) + abs(estimated_haplotype_matrix[marker*total_haplotypes+current_subject*2+1]-h2);
  //cerr<<"Diff is "<<diff<<endl;
  //diff = abs(h1+h2-estimated_haplotype_matrix[marker*total_haplotypes+current_subject*2]-estimated_haplotype_matrix[marker*total_haplotypes+current_subject*2+1]);
  //cerr<<"Diff2 is "<<diff<<endl;
  site_matches[marker] += 2-diff;
  site_mismatches[marker]+=diff;
  //cerr<<"Imputed haplotype at marker: "<<marker<<" from "<<hap1<<","<<hap2<<endl;
}

void HMM_impute::compute_genotype_probs(int marker_loc, int & offset){
}

void HMM_impute::extract_backward_prob(int marker_loc, int & offset, int hap1, int hap2, float *  marginals){
  double  start = clock();
  int block_id = (marker_loc / matrix_compression_ratio);
  offset = block_id * matrix_compression_ratio;
  int block_marker_end;
  int signpost_row = (matrix_signposts-1)-block_id;
  if (block_id!=current_block_id){
    //cerr<<"Opening inputfile\n";
    //cerr<<"Backward rescalings in extract block: "<<backward_rescalings<<endl;
    ifs_cache.open(cachefilename.data(),ios::in|ios::binary);
    ofs<<"Opened inputfile\n";
    current_block_id = block_id;
    int marker_start = block_id*matrix_compression_ratio;
    int marker_end = marker_start+matrix_compression_ratio; 
    if (marker_end>markers) marker_end = markers;
    block_marker_end = marker_end - marker_start;
    if(settings->use_gpu){
      #ifdef USE_GPU
      err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.marker_offset_buffer, CL_TRUE, 0,  sizeof(int)*1, &offset , NULL, NULL );
      clSafe(err, "write marker offset index");
      ofs<<"in extract backward prob, fetching signpost row "<<signpost_row<<endl;
      if(use_compressed){
        flush2ram(ifs_cache,short_backward_mat_block + (block_marker_end-1)*compressed_states,compressed_states,signpost_row);
        err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.backward_block_buffer, CL_TRUE, (block_marker_end-1)*compressed_states*sizeof(packedshort_t),  compressed_states*sizeof(packedshort_t), short_backward_mat_block+(block_marker_end-1)*compressed_states, NULL, NULL );
      }else{
        flush2ram(ifs_cache,backward_mat_block + (block_marker_end-1)*states,states,signpost_row);
        err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.backward_block_buffer, CL_TRUE, (block_marker_end-1)*states*sizeof(float), states*sizeof(float), backward_mat_block+(block_marker_end-1)*states, NULL, NULL );
      }
      clSafe(err, "write backward block");
      start = clock();
      int block_marker;
      int marker;
      block_marker = block_marker_end-1;
      marker = offset+block_marker;
      if (use_newfb){
        err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.marker_index_buffer, CL_TRUE, 0,  sizeof(int)*1, &marker, NULL, NULL );
        clSafe(err, "write backward marker index");
        init_forwardbackward_gpu();
      }
      bool debug_rescalings = false;
      for(block_marker = block_marker_end-2;block_marker>=0;--block_marker){
        marker = offset+block_marker;
        //cerr<<"Marker: "<<marker<<endl;
        err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.marker_index_buffer, CL_TRUE, 0,  sizeof(int)*1, &marker, NULL, NULL );
        clSafe(err, "write backward marker");
        forwardbackward_gpu(marker);
        if(debug_rescalings){
          err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.backward_underflows_buffer, CL_TRUE, 0, 1*sizeof(int),&backward_rescalings);
          clSafe(err, "read rescalings");
          ofs<<"extract block marker "<<marker<<", backward rescalings: "<<backward_rescalings<<endl;
        }
      }
      ofs<<"compute backward block time: "<<(clock()-start)/CLOCKS_PER_SEC<<endl;
      if (block_id==0){
        // compute log likelihood
        if(use_compressed){
          err = opencl_info.command_queue.enqueueNDRangeKernel(opencl_info.kernel_loglikelihood_short,cl::NullRange,cl::NDRange(BLOCK_WIDTH,1),cl::NDRange(BLOCK_WIDTH,1),NULL,NULL);
        }else{
          err = opencl_info.command_queue.enqueueNDRangeKernel(opencl_info.kernel_loglikelihood_float,cl::NullRange,cl::NDRange(BLOCK_WIDTH,1),cl::NDRange(BLOCK_WIDTH,1),NULL,NULL);
        }
        clSafe(err,"launch loglikelihood");
        err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.backward_underflows_buffer, CL_TRUE, 0, 1*sizeof(int),&backward_rescalings);
        clSafe(err, "read rescalings");
        ofs<<"backward rescalings: "<<backward_rescalings<<endl;
        err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.loglikelihood_buffer, CL_TRUE, 0, 1*sizeof(float),&likelihood_log);
        clSafe(err, "read log likelihood");
        ofs<<"Log likelihood: "<<likelihood_log<<endl;
        cerr<<"Log likelihood: "<<likelihood_log<<endl;
      }
      #endif
    }else{
      // seed the backward probability vector
      //ofs<<"extract block, marker "<<offset+block_marker_end-1<<", fetching signpost row "<<signpost_row<<endl;
      if(use_compressed){
        flush2ram(ifs_cache,short_backward_mat_block + (block_marker_end-1)*compressed_states,compressed_states,signpost_row);
      }else{
        flush2ram(ifs_cache,backward_mat_block + (block_marker_end-1)*states,states,signpost_row);
      }
      // Compute backward probabilities
      double start = clock();
      for(int block_marker = block_marker_end-2;block_marker>=0;--block_marker){
        int marker = offset+block_marker;
        //ofs<<"extract block backward marker "<<marker<<endl;
        if(use_compressed){
          forwardbackward_cpu(marker,short_backward_mat_block+(block_marker+1)*
          compressed_states,short_backward_mat_block+(block_marker)*
          compressed_states,site_thetas[marker+1],backward_rescalings);
        }else{
          forwardbackward_cpu(marker, backward_mat_block+(block_marker+1)*
          states,backward_mat_block+(block_marker)*states,
          site_thetas[marker+1],backward_rescalings);
        }
      }
      if (block_id==0){
        if(use_compressed){
          short2floatarr(short_backward_mat_block,temp_prob_vec);
        } 
        float normalizer = 0;
        for(int hap1=0;hap1<total_template_haplotypes;++hap1){
          int state_offset = hap1*total_template_haplotypes;
          for(int hap2=0;hap2<=hap1;++hap2){
            int state = state_offset+hap2;
        // the following contains the last marker from backward algorithm and the first marker from the forward algorithm
            if(use_compressed){
              normalizer += temp_prob_vec[state] * prior_vec[state];
            }else{
              normalizer += backward_mat_block[state] * prior_vec[state];
            }
          }
        } 
        likelihood_log = log(normalizer) - (backward_rescalings)*log(LARGE);
        cerr<<"backward rescalings: "<<backward_rescalings<<endl;
        cerr<<"Log likelihood: "<<likelihood_log<<endl;
      }
      cerr<<"Extract block time: "<<(clock()-start)/CLOCKS_PER_SEC<<endl; 
      // end of backward prob 
    }
    //cerr<<"Closing input cache\n";
    ifs_cache.close();
    ofs<<"Closed input cache\n";
    // now compute the marginals
  }
  if(settings->use_gpu){
    #ifdef USE_GPU
    int hap_pair[] = {hap1,hap2};
    int block_marker = (marker_loc-offset);
    err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.hap_pair_buffer, CL_TRUE, 0,  sizeof(int)*2, hap_pair , NULL, NULL );
    clSafe(err, "write hap_pair");
    err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.marker_index_buffer, CL_TRUE, 0,  sizeof(int)*1, &block_marker , NULL, NULL );
    clSafe(err, "write block row");
    cl::Event * event ;
    bool benchmark = false;
    if (benchmark){
      event = new cl::Event;
    }else{
      event = NULL;
    }
//memset(marginals,1,sizeof(float)*4);
//return;
    double start = clock();
    if(use_compressed){
      err = opencl_info.command_queue.enqueueNDRangeKernel(opencl_info.kernel_collect_xover_marginals_block_short,cl::NullRange,cl::NDRange(BLOCK_WIDTH*total_template_haplotypes,1),cl::NDRange(BLOCK_WIDTH,1),NULL,event);
    }else{
      err = opencl_info.command_queue.enqueueNDRangeKernel(opencl_info.kernel_collect_xover_marginals_block_float,cl::NullRange,cl::NDRange(BLOCK_WIDTH*total_template_haplotypes,1),cl::NDRange(BLOCK_WIDTH,1),NULL,event);
    }
    clSafe(err,"launch xover marginals");
    err = opencl_info.command_queue.enqueueNDRangeKernel(opencl_info.kernel_collect_xover_marginals2_block,cl::NullRange,cl::NDRange(BLOCK_WIDTH,4),cl::NDRange(BLOCK_WIDTH,1),NULL,event);
    clSafe(err,"launch xover marginals");
    if (benchmark){
      event->wait();
      cerr<<"Marker: "<<marker_loc<<": Collect xover marginals: "<<(clock()-start)<<endl;
      delete event;
    }
    err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.xover_full_marginal_buffer, CL_TRUE, 0, (4)*sizeof(float),marginals);
    clSafe(err, "read xover marginals");
    bool debug_xover_marginals = false;
    if (debug_xover_marginals){
      int block_marker = (marker_loc-offset);

      float big_marginals[4*total_template_haplotypes];
      err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.xover_haplotype_marginals_buffer, CL_TRUE, 0, (4*total_template_haplotypes)*sizeof(float),big_marginals);
      clSafe(err, "read xover big marginals");
      ofs<<"gpu marginals for marker: "<<marker_loc<<"\n";
      for(int i=0;i<4;++i){
        ofs<<i<<":"<<marginals[i]<<endl;
      }
      ofs<<"retrieving backward block row\n";
      float * blockptr;
      if(use_compressed){
        err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.backward_block_buffer, CL_TRUE, block_marker*compressed_states*sizeof(packedshort_t), compressed_states*sizeof(packedshort_t),short_backward_mat_block + block_marker*compressed_states);
        short2floatarr(short_backward_mat_block+(block_marker)*
        compressed_states, temp_prob_vec);
        blockptr = temp_prob_vec;
      }else{
        err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.backward_block_buffer, CL_TRUE, block_marker*states*sizeof(float), states*sizeof(float),backward_mat_block + block_marker*states);
        blockptr = backward_mat_block+block_marker*states;
      }
      clSafe(err, "read backwardvec");
      float marginals2[4];
      memset(marginals2,0,sizeof(float)*4);
      for(int state=0;state<states;++state){
        if (state_haplo_indices[2*state]<=state_haplo_indices[2*state+1]){
          //blockptr[state] = 1;
          // for float recombinants any haplotype qualifies
          marginals2[MARGINALS_BOTH_XOVER] += blockptr[state];
          // for first recombinant
          if(hap2==state_haplo_indices[state*2] ||
          hap2==state_haplo_indices[state*2+1])
            marginals2[MARGINALS_FIRST_XOVER] += blockptr[state];
          // for second recombinant
          if(hap1==state_haplo_indices[state*2] ||
          hap1==state_haplo_indices[state*2+1])
            marginals2[MARGINALS_SECOND_XOVER] += blockptr[state];
          // for no crossover
          if(hap1==state_haplo_indices[state*2] &&
          hap2==state_haplo_indices[state*2+1] ||
          hap2==state_haplo_indices[state*2] &&
          hap1==state_haplo_indices[state*2+1] )
            marginals2[MARGINALS_NO_XOVER] += blockptr[state];
        }
      }
      ofs<<"cpu marginals for marker: "<<marker_loc<<"\n";
      for(int i=0;i<4;++i){
        ofs<<marginals2[i]<<endl;
      }
      cerr<<"Wrote xover marginals\n";
      MPI_Finalize();
      //exit(0);
    }
    #endif
  }else{
    // Compute the 4 marginals
    //cerr<<"Getting blockrow: "<<marker_loc-offset<<endl;
    //float * backward_vec = backward_mat_block+(marker_loc-offset)*states; 
    //ofs<<"Computing marginals for marker "<<marker_loc<<endl;
    float * blockptr;
    if(use_compressed){
      short2floatarr(short_backward_mat_block+(marker_loc-offset)*
      compressed_states, temp_prob_vec);
      blockptr = temp_prob_vec;
    }else{
      blockptr = backward_mat_block+(marker_loc-offset)*states;
    }
    memset(marginals,0,sizeof(float)*4);
    for(int hapindex1=0;hapindex1<total_template_haplotypes;++hapindex1){
      for(int hapindex2=0;hapindex2<=hapindex1;++hapindex2){
        int state = hapindex1*total_template_haplotypes+hapindex2;
        // for float recombinants any haplotype qualifies
        marginals[MARGINALS_BOTH_XOVER] += blockptr[state];
        // for first recombinant
        if(hap2==hapindex1 || hap2==hapindex2){ 
          marginals[MARGINALS_FIRST_XOVER] += blockptr[state];
        }
        // for second recombinant
        if(hap1==hapindex1 || hap1==hapindex2) 
          marginals[MARGINALS_SECOND_XOVER] += blockptr[state];
        // for no crossover
        if((hap1==hapindex1 && hap2==hapindex2) || 
        (hap2==hapindex1 && hap1==hapindex2)) 
          marginals[MARGINALS_NO_XOVER] += blockptr[state];
      }
    }
  }
  for(int i=0;i<4;++i){
     //cerr<<marginals[i]<<endl;
  }
  //exit(0);
  return;
}

void HMM_impute::compute_genotype_posterior(int marker,int offset){
}

int HMM_impute::randomstate(int marginal_type,int block_row,int hap1,
int hap2,float random){
  int state = 0;
  if (trace)ofs<<"Requesting marginaltype: "<<marginal_type<<", blockrow: "<<block_row<<", haps: "<<hap1<<","<<hap2<<", random: "<<random<<endl;
  float * blockptr;
  if (marginal_type!=MARGINALS_NO_XOVER){
    if (settings->use_gpu){
#ifdef USE_GPU
      if(use_compressed){
        err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.backward_block_buffer, CL_TRUE, block_row*compressed_states*sizeof(packedshort_t), compressed_states*sizeof(packedshort_t),short_backward_mat_block + block_row*compressed_states,NULL,NULL);
        clSafe(err, "read backwardvec");
      }else{
        err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.backward_block_buffer, CL_TRUE, block_row*states*sizeof(float), states*sizeof(float),backward_mat_block + block_row*states,NULL,NULL);
        clSafe(err, "read backwardvec");
      }
      //ostringstream oss1;
      //oss1<<"debug.backward.block."<<block_row<<".mpirank."<<mpi_rank;
      //debugfloat(oss1.str().data(),backward_mat_block + block_row*states,states);
#endif
    }
     
    if(use_compressed){ 
      short2floatarr(short_backward_mat_block+block_row*compressed_states,
      temp_prob_vec);
      blockptr = temp_prob_vec;
    }else{
      blockptr = backward_mat_block+block_row*states;
    }
  }
  float sum = 0;
  bool found = false;
  //int hapindex1=0,hapindex2=0;
  float * hap1ptr = blockptr+hap1*total_template_haplotypes;
  float * hap2ptr = blockptr+hap2*total_template_haplotypes;
  int hapiter = 0;
  //marginal_type = MARGINALS_NO_XOVER;
  switch(marginal_type){
    case MARGINALS_NO_XOVER:
      state = (hap1>hap2)?hap1*total_template_haplotypes+hap2:
      hap2*total_template_haplotypes+hap1;
      found = true;
      break;
    case MARGINALS_FIRST_XOVER:
      //for(int o = hap2+1;o<total_template_haplotypes;++o){
      //  hap2ptr[o] = temp_prob_vec[o*total_template_haplotypes+hap2];
      //}
      while(!found && hapiter<total_template_haplotypes){
        sum+=hap2ptr[hapiter++];
        found = random<sum;  
      }
      --hapiter;
      state = (hap2<hapiter)?hapiter*total_template_haplotypes+hap2:
      hap2*total_template_haplotypes+hapiter;
      break;
    case MARGINALS_SECOND_XOVER:
      //for(int o = hap1+1;o<total_template_haplotypes;++o){
      //  hap1ptr[o] = temp_prob_vec[o*total_template_haplotypes+hap1];
      //}
      while(!found && hapiter<total_template_haplotypes){
        sum+=hap1ptr[hapiter++];
        found = random<sum;  
      }
      --hapiter;
      state = (hap1<hapiter)?hapiter*total_template_haplotypes+hap1:
      hap1*total_template_haplotypes+hapiter;
      break;
    case MARGINALS_BOTH_XOVER:
      for(int hap1=0;hap1<total_template_haplotypes;++hap1){
        int state_offset = hap1*total_template_haplotypes;
        for(int hap2=0;hap2<=hap1;++hap2){
          state = state_offset+hap2;
          sum+=blockptr[state];
          found = random<sum;
          if (found) break;
        }
        if (found) break;
      }
      break;
  }
  last_block_row_retrieved = block_row;
  if (!found){ 
    cerr<<"Marginal type is "<<marginal_type<<endl;
    cerr<<"Failed to find a valid state. Random/sum: "<<random<<","<<sum<<" Retry.\n";
    state = states;
  }
  return state;
}

void HMM_impute::sample_chromosomes(){
  // choose two template haplotypes at random
  last_block_row_retrieved = -1;
  int offset;
  float marginals[4];
  extract_backward_prob(0,offset,0,0,marginals);
  for(int i=0;i<4;++i){
    //cerr<<"Marginals: "<<marginals[i]<<endl;
  }
  float normalizer = marginals[MARGINALS_BOTH_XOVER];
  int state = 0;
  do{
    float r = math->RandomUniform()*normalizer;
    state = randomstate(MARGINALS_BOTH_XOVER,0,0,0,r);
    //cerr<<"Normalizer: "<<normalizer<<", r: "<<r<<" ,state: "<<state<<endl;
  }while(state>=states);
  int hap1 = state_haplo_indices[state*2];
  int hap2 = state_haplo_indices[state*2+1];
  ofs<<"Init Choosing haplotypes: "<<hap1<<","<<hap2<<endl;
  impute_and_evaluate(0,hap1,hap2);
  // from these two random templates, copy them to the estimated_haplotype
  for(int marker=1;marker<markers;++marker){
    // for interval between j and j+1 determine how many crossovers to insert
    // we will have 00,01,10,11 as possibilities
    extract_backward_prob(marker,offset,hap1,hap2,marginals);
    for(int i=0;i<4;++i){
      //cerr<<"Marginals "<<i<<": "<<marginals[i]<<endl;
    }
    float theta = site_thetas[marker];
    float no_xover = marginals[MARGINALS_NO_XOVER] * pow(1.-theta,2);
    float xover1 = marginals[MARGINALS_FIRST_XOVER] * theta * (1-theta) / 
    total_haplotypes;
    float xover2 = marginals[MARGINALS_SECOND_XOVER] * theta * (1-theta) / 
    total_haplotypes;
    float two_xover = marginals[MARGINALS_BOTH_XOVER] * pow(theta,2) / 
    (total_haplotypes*total_haplotypes);
    float r = math->RandomUniform()*(no_xover+xover1+xover2+two_xover);
    if (r < no_xover){  // no crossover
      // copy both template haplotypes to the estimated haplotype
      // do nothing  
      //cerr<<"No crossover"<<endl;
    }else if (r < no_xover+xover1){
      // choose a new value for hap1 and fill in
      do{
        float r = math->RandomUniform() * marginals[MARGINALS_FIRST_XOVER];
        state = randomstate(MARGINALS_FIRST_XOVER,marker-offset,hap1,hap2,r);
        if (state==states){
          ofs<<"repeat firstxover.";
          for(int i=0;i<4;++i){
            ofs<<" "<<marginals[i];
          }
          ofs<<endl;
        }
      }while(state>=states);
      ++site_xovers[marker];
      hap1 = state_haplo_indices[state*2]==hap1?state_haplo_indices[state*2+1]:state_haplo_indices[state*2];
      ofs<<"1st crossover"<<endl;
      ofs<<"Choosing haplotypes: "<<hap1<<","<<hap2<<endl;
    }else if (r<no_xover+xover1+xover2){
      // choose a new value for hap2 and fill in
      do{
        float r = math->RandomUniform() * marginals[MARGINALS_SECOND_XOVER];
        state = randomstate(MARGINALS_SECOND_XOVER,marker-offset,hap1,hap2,r);
        if (state==states){
          ofs<<"repeat seconddxover.";
          for(int i=0;i<4;++i){
            ofs<<" "<<marginals[i];
          }
          ofs<<endl;
        }
      }while(state>=states);
      ++site_xovers[marker];
      hap2 = state_haplo_indices[state*2+1]==hap2?state_haplo_indices[state*2]:state_haplo_indices[state*2+1];
      ofs<<"2nd crossover"<<endl;
      ofs<<"Choosing haplotypes: "<<hap1<<","<<hap2<<endl;
    }else{
      // choose a new value for hap1 and fill in
      // choose a new value for hap2 and fill in
      do{
        float r = math->RandomUniform() * marginals[MARGINALS_BOTH_XOVER];
        state = randomstate(MARGINALS_BOTH_XOVER,marker-offset,hap1,hap2,r);
        if (state==states){
          ofs<<"repeat bothxover.";
          for(int i=0;i<4;++i){
            ofs<<" "<<marginals[i];
          }
          ofs<<endl;
        }
      }while(state>=states);
      hap1 = state_haplo_indices[state*2];
      hap2 = state_haplo_indices[state*2+1];
      ofs<<"Both crossover"<<endl;
      ofs<<"Choosing haplotypes: "<<hap1<<","<<hap2<<endl;
      site_xovers[marker]+=2;
    }
    //cerr<<"At marker: "<<marker<<", haplo indices are: "<<hap1<<","<<hap2<<endl;
    impute_and_evaluate(marker,hap1,hap2);
  }
  //cerr<<"Sample chr time: "<<(clock()-start)/CLOCKS_PER_SEC<<endl;
}

void HMM_impute::forwardbackward_cpu(int marker,packedshort_t * short_previous_prob_vec,packedshort_t * short_current_prob_vec, float theta, int & rescaling){
  // EMULATE THE COMPUTE_EMISSION KERNEL
  ofs<<"Emulating emission compute and marginals kernel\n";
  short2floatarr(short_previous_prob_vec,temp_prob_vec);
  compute_emission_vec(marker+1,emission_vec);
  bool debug_emission = true;
  if (debug_emission){
    ostringstream oss;
    cerr<<"Reading from CPU at marker "<<marker+1<<endl;
    oss<<"debug.cpu.emission.marker."<<marker+1<<".mpirank."<<mpi_rank;
    debugfloat(oss.str().data(),emission_vec,states);
    exit(0);
  }
  for(int state=0;state<states;++state){
    temp_prob_vec[state] *= emission_vec[state];
  }
  // compute the marginals
  memset(haplotype_marginals,0,sizeof(float)*2*total_template_haplotypes);
  bool odd = marker%2;
  float joint_marginal2=0;
  for(int hap1=0;hap1<total_template_haplotypes;++hap1){
    float offdiag_sum=0;
    int state_offset = hap1*total_template_haplotypes;
    for(int hap2=0;hap2<hap1;++hap2){
      float prev = temp_prob_vec[state_offset+hap2];
      haplotype_marginals[odd*total_template_haplotypes+hap2]+=prev;
      offdiag_sum+=prev;
    }
    float other = temp_prob_vec[state_offset+hap1]+offdiag_sum;
    haplotype_marginals[odd*total_template_haplotypes+hap1]+=other;
    joint_marginal2+=other+offdiag_sum;
  }
  bool debug_marginals = true;
  if (debug_marginals){
    ostringstream oss;
    oss<<"debug.cpu.marginals.mpirank."<<mpi_rank;
    debugfloat(oss.str().data(),haplotype_marginals+odd*total_template_haplotypes,total_template_haplotypes);
    cerr<<"Joint marginal: "<<joint_marginal2<<endl;
    MPI_Finalize();
  }
  // this represents a store into global memory
  float2shortarr(temp_prob_vec,short_temp_prob_vec);
  // EMULATE THE FORWARD BACKWARD KERNEL
  ofs<<"Emulating forward backward kernel\n";
  // this represents a load from global memory
  short2floatarr(short_temp_prob_vec,temp_prob_vec);
  float change0 = (1-theta)*(1-theta);
  float change1 = (1-theta)*theta/(total_template_haplotypes);
  float change2 = joint_marginal2*(theta)*(theta)/(total_template_haplotypes*total_template_haplotypes);
  //if (joint_marginal2<1e-2){
  if (joint_marginal2<SMALL){
    change0*=LARGE;
    change1*=LARGE;
    change2*=LARGE;
    //cerr<<"Rescaled at marker "<<marker<<"!\n";
    ++rescaling;
  }
  for(int haploindex=0;haploindex<total_template_haplotypes;++haploindex){
    int state_offset = haploindex*total_template_haplotypes;
    for(int haploindex2=0;haploindex2<haploindex; ++haploindex2){
      int state = state_offset+haploindex2;
      temp_prob_vec[state] = temp_prob_vec[state]*change0+change1*(haplotype_marginals[odd*total_template_haplotypes+haploindex]+haplotype_marginals[odd*total_template_haplotypes+haploindex2])+2*change2;
      int otherstate = haploindex2*total_template_haplotypes+haploindex;
      temp_prob_vec[otherstate] = temp_prob_vec[state] ;
    }
    temp_prob_vec[state_offset+haploindex] = temp_prob_vec[state_offset+haploindex]*change0+change1 * haplotype_marginals[odd*total_template_haplotypes+haploindex] + change2;
  }
  bool debug_fb = false;
  if (debug_fb) {
    ostringstream oss; 
    oss<<"debug.cpu.fb."<<marker<<".mpirank."<<mpi_rank;
    debugfloat(oss.str().data(),temp_prob_vec,states);
    cerr<<"Wrote FB to disk\n";
    MPI_Finalize();
  }
  //store in global memory
  float2shortarr(temp_prob_vec,short_current_prob_vec);
}

void HMM_impute::forwardbackward_cpu(int marker,float * previous_prob_vec,float * current_prob_vec, float theta, int & rescaling){
  //ofs<<"Marker "<<marker<<endl;
//  bool debug_prior = false;
//        if (debug_prior) {
//          ostringstream oss; 
//          oss<<"debug.signpost.prior."<<marker<<".mpirank."<<mpi_rank;
//          debugfloat(oss.str().data(),previous_prob_vec,states);
//          exit(0);
//        }
//      
  compute_emission_vec(marker+1,emission_vec);
  bool debug_emission = false;
  if (debug_emission){
    ostringstream oss;
    cerr<<"Reading from CPU at marker "<<marker+1<<endl;
    oss<<"debug.cpu.emission.marker."<<marker+1<<".mpirank."<<mpi_rank;
    debugfloat(oss.str().data(),emission_vec,states);
    exit(0);
  }
  for(int i=0;i<states;++i){
    temp_prob_vec[i] = previous_prob_vec[i] * emission_vec[i];
  }
//  // initialize marginal emissions to zero
  //cerr<<"clearing hap marginals\n";
  memset(haplotype_marginals,0,sizeof(float)*2*total_template_haplotypes);
  //cerr<<"cleared hap marginals\n";
  float full_marginal=0;
  bool odd = marker%2;
  for(int hap1=0;hap1<total_template_haplotypes;++hap1){
    float offdiag_sum=0;
    int state_offset = hap1*total_template_haplotypes;
    for(int hap2=0;hap2<hap1;++hap2){
      float prev = temp_prob_vec[state_offset+hap2];
      haplotype_marginals[odd*total_template_haplotypes+hap2]+=prev;
      offdiag_sum+=prev;
    }
    float other = temp_prob_vec[state_offset+hap1]+offdiag_sum;
    haplotype_marginals[odd*total_template_haplotypes+hap1]+=other;
    full_marginal+=other+offdiag_sum;
  }
  bool debug_marginals = false;
  if (debug_marginals){
    ostringstream oss;

    oss<<"debug.cpu.marginals.mpirank."<<mpi_rank;
    debugfloat(oss.str().data(),haplotype_marginals+odd*total_template_haplotypes,total_template_haplotypes);
    cerr<<"Full marginal: "<<full_marginal<<endl;
    MPI_Finalize();
    //exit(0);
  }
  float change0 = (1-theta)*(1-theta);
  float change1 = (1-theta)*theta/(total_template_haplotypes);
  float change2 = full_marginal*(theta)*(theta)/(total_template_haplotypes*total_template_haplotypes);
  if (full_marginal<SMALL){
    change0*=LARGE;
    change1*=LARGE;
    change2*=LARGE;
    //cerr<<"Rescaled at marker "<<marker<<"!\n";
    ++rescaling;
  }
  for(int haploindex=0;haploindex<total_template_haplotypes;++haploindex){
    int state_offset = haploindex*total_template_haplotypes;
    for(int haploindex2=0;haploindex2<haploindex; ++haploindex2){
      int state = state_offset+haploindex2;
      current_prob_vec[state] = temp_prob_vec[state]*change0+change1*(haplotype_marginals[odd*total_template_haplotypes+haploindex]+haplotype_marginals[odd*total_template_haplotypes+haploindex2])+2*change2;
      int otherstate = haploindex2*total_template_haplotypes+haploindex;
      current_prob_vec[otherstate] = current_prob_vec[state] ;
    }
    current_prob_vec[state_offset+haploindex] = temp_prob_vec[state_offset+haploindex]*change0+change1 * haplotype_marginals[odd*total_template_haplotypes+haploindex] + change2;
  }
  bool debug_fb = marker==-7986;
        if (debug_fb) {
//          float minlog = 100,maxlog = -100;
//          for(int i=0;i<states;++i){
//            float lprob =  round(log2(current_prob_vec[i]));
//            if (lprob>maxlog) maxlog = lprob;
//            if (lprob<minlog) minlog = lprob;
//          }
//          ofs<<"marker: "<<marker<<", min: "<<minlog<<", maxlog: "<<maxlog<<endl;
          ostringstream oss; 
          oss<<"debug.cpu.fb."<<marker<<".mpirank."<<mpi_rank;
//          ofstream ofs1(oss.str().data());
//          for(int i=0;i<1000;++i){
//            float newval = short2float(float2short(current_prob_vec[i]));
//            ofs1<<i<<"\t"<<current_prob_vec[i]<<"\t"<<newval<<endl;
//          }
//          ofs1.close();
          if(use_newfb){
            compute_emission_vec(marker,emission_vec);
            for(int i=0;i<states;++i){
              temp_prob_vec[i] = current_prob_vec[i] * emission_vec[i];
            }
            debugfloat(oss.str().data(),temp_prob_vec,states);
          }else{
            debugfloat(oss.str().data(),current_prob_vec,states);
          }
          cerr<<"Wrote CPU fb debug\n";
          MPI_Finalize();
//          //exit(0);
        }
  return;
}

void HMM_impute::init_forwardbackward_gpu(){
  bool benchmark = false;
  cl::Event * event;
  if (benchmark){
    event = new cl::Event;
  }else{
    event = NULL;
  }
  double start;
  start = clock();
  if (use_compressed){
    err = opencl_info.command_queue.enqueueNDRangeKernel(opencl_info.kernel_init_forwardbackward_short,cl::NullRange,cl::NDRange(total_template_haplotypes*BLOCK_WIDTH,1),cl::NDRange(BLOCK_WIDTH,1),NULL,event);
  }else{ 
    err = opencl_info.command_queue.enqueueNDRangeKernel(opencl_info.kernel_init_forwardbackward_float,cl::NullRange,cl::NDRange(total_template_haplotypes*BLOCK_WIDTH,1),cl::NDRange(BLOCK_WIDTH,1),NULL,event);
  }
  clSafe(err,"launch init_forwardbackward");
  if(benchmark){
    event->wait();
    cerr<<"emission vec: "<<clock()-start<<endl;
  }
  bool debug_emission = false;
  if (debug_emission){
    err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.emission_vec_buffer, CL_TRUE, 0, (states)*sizeof(float),emission_vec);
    clSafe(err, "read emission_vec backward");
    ostringstream oss; 
    oss<<"debug.gpu.emission.first.mpirank."<<mpi_rank;
    debugfloat(oss.str().data(),emission_vec,states);
    cerr<<"Read from GPU\n";
    MPI_Finalize();
  }
}

void HMM_impute::forwardbackward_gpu(int marker){
  bool benchmark = false;
  cl::Event * event;
  if (benchmark){
    event = new cl::Event;
  }else{
    event = NULL;
  }
  double start;
  if (!use_newfb){
    start = clock();
    if (use_compressed){
      err = opencl_info.command_queue.enqueueNDRangeKernel(opencl_info.kernel_init_forwardbackward_short,cl::NullRange,cl::NDRange(total_template_haplotypes*BLOCK_WIDTH,1),cl::NDRange(BLOCK_WIDTH,1),NULL,event);
    }else{ 
      err = opencl_info.command_queue.enqueueNDRangeKernel(opencl_info.kernel_init_forwardbackward_float,cl::NullRange,cl::NDRange(total_template_haplotypes*BLOCK_WIDTH,1),cl::NDRange(BLOCK_WIDTH,1),NULL,event);
    }
    clSafe(err,"launch compute_emission_vector");
    if(benchmark){
      event->wait();
      cerr<<"emission vec: "<<clock()-start<<endl;
    }
    bool debug_emission = false;
    if (debug_emission){
      err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.emission_vec_buffer, CL_TRUE, 0, (states)*sizeof(float),emission_vec);
      cerr<<"mpirankA:"<<mpi_rank<<endl;
      clSafe(err, "read emission_vec backward");
      cerr<<"mpirankB:"<<mpi_rank<<endl;
      
      ostringstream oss; 
      oss<<"debug.gpu.emission.first.mpirank."<<mpi_rank;
      debugfloat(oss.str().data(),emission_vec,states);
      bool odd = use_newfb?marker%2:0;
      err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.haplotype_marginals_buffer, CL_TRUE, odd*total_template_haplotypes*sizeof(float), total_template_haplotypes*sizeof(float),haplotype_marginals);
      clSafe(err, "read backward marginals");
      ostringstream oss1;
      oss1<<"debug.gpu.marginals.mpirank."<<mpi_rank;
      debugfloat(oss1.str().data(),haplotype_marginals,total_template_haplotypes);
      cerr<<"Read from GPU at marker "<<marker+1<<"\n";
      MPI_Finalize();
    }
  }
  start = clock();
  err = opencl_info.command_queue.enqueueNDRangeKernel(opencl_info.kernel_reduce_fb_marginals,cl::NullRange,cl::NDRange(BLOCK_WIDTH),cl::NDRange(BLOCK_WIDTH),NULL,event);
  clSafe(err,"launch reduce_fb_marginals");
  if(benchmark){
    event->wait();
    cerr<<"Reduce marginals: "<<clock()-start<<endl;
  }
  bool debug_marginals = false;
  if (debug_marginals){
    float full_marginal;
    err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.full_marginal_buffer, CL_TRUE, 0, sizeof(float),&full_marginal);
    clSafe(err, "read backward marginal");
    cerr<<"Joint marginal: "<<full_marginal<<endl;
    MPI_Finalize();
  }
  start = clock();
  if(use_compressed){
    err = opencl_info.command_queue.enqueueNDRangeKernel(opencl_info.kernel_forward_backward_short,cl::NullRange,cl::NDRange(total_template_haplotypes*BLOCK_WIDTH,1),cl::NDRange(BLOCK_WIDTH,1),NULL,event);
  }else{
    if (use_newfb){
      err = opencl_info.command_queue.enqueueNDRangeKernel(opencl_info.kernel_forward_backward_emission_float,cl::NullRange,cl::NDRange(total_template_haplotypes*BLOCK_WIDTH,1),cl::NDRange(BLOCK_WIDTH,1),NULL,event);
    }else{
      err = opencl_info.command_queue.enqueueNDRangeKernel(opencl_info.kernel_forward_backward_float,cl::NullRange,cl::NDRange(total_template_haplotypes*BLOCK_WIDTH,1),cl::NDRange(BLOCK_WIDTH,1),NULL,event);
    }
  }
  clSafe(err,"launch forwardbackward");
  if(benchmark){
    event->wait();
    cerr<<"forward backward: "<<clock()-start<<endl;
  }

  bool debug_fb = marker==-7988;
  if (debug_fb){
    ostringstream oss1; 
    if (use_compressed){ 
      err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.emission_vec_buffer, CL_TRUE, 0, (states)*sizeof(float),emission_vec);
      clSafe(err, "read emission_vec backward");
      oss1<<"debug.gpu.fb."<<marker<<".mpirank."<<mpi_rank;
      debugfloat(oss1.str().data(),emission_vec,states);
    }else{
      err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.current_prob_vec_buffer, CL_TRUE, 0, (states)*sizeof(float),current_prob_vec);
      clSafe(err, "read current prob vec backward");
      oss1<<"debug.gpu.fb."<<marker<<".mpirank."<<mpi_rank;
      debugfloat(oss1.str().data(),current_prob_vec,states);
    }
    cerr<<"Wrote FB to disk\n";
    MPI_Finalize();
  }
  if (benchmark) delete event;
  return;
}



void HMM_impute::compute_matrix_signposts(){
  int signpost = 0;// set the sign post to the last one left off by the forward alg
  //cerr<<"Signpost initialized: "<<signpost<<endl;
  ofs_cache.open(cachefilename.data(),ios::out|ios::binary);
  ofs<<"Opened outputfile\n";
  if(use_compressed){
    flush2disk(ofs_cache,short_prior_vec,compressed_states);
  }else{
    flush2disk(ofs_cache,prior_vec,states);
  }
  ++signpost;
  backward_rescalings = 0;
  if (settings->use_gpu){
    #ifdef USE_GPU
    int offset = -1;
    err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.marker_offset_buffer, CL_TRUE, 0,  sizeof(int)*1, &offset , NULL, NULL );
    clSafe(err, "write marker offset index");
    // write the last backward prob vec
    if(use_compressed){
      err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.current_prob_vec_buffer, CL_TRUE, 0,  sizeof(packedshort_t)*compressed_states, short_prior_vec , NULL, NULL );
    }else{
      err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.current_prob_vec_buffer, CL_TRUE, 0,  sizeof(float)*states, prior_vec , NULL, NULL );
    }
    clSafe(err, "write prior vec");
    // store the start and last signpost
    err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.backward_underflows_buffer, CL_TRUE, 0,  sizeof(int)*(1), &backward_rescalings , NULL, NULL );
    clSafe(err, "write buffer backward under start");
    int backward_marker;
    int  writesignpost = 0;
    if (use_newfb){
      backward_marker = markers - 1;
      err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.marker_index_buffer, CL_TRUE, 0,  sizeof(int)*1, &backward_marker , NULL, NULL );
      clSafe(err, "write backward marker index");
      init_forwardbackward_gpu();
      err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.writesignpost_buffer, CL_TRUE, 0,  sizeof(int)*1, &writesignpost , NULL, NULL );
      clSafe(err, "setting signpost write to 0");
    }
    bool debug_rescalings = false;
    for(backward_marker = markers-2;backward_marker>=(matrix_compression_ratio-1);--backward_marker){
      //ofs<<"Compute signpost marker: "<<backward_marker<<endl;
      err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.marker_index_buffer, CL_TRUE, 0,  sizeof(int)*1, &backward_marker , NULL, NULL );
      clSafe(err, "write backward marker index");
      writesignpost = (backward_marker+1) % matrix_compression_ratio == 0;
      if(use_newfb){
        if (writesignpost){
          err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.writesignpost_buffer, CL_TRUE, 0,  sizeof(int)*1, &writesignpost , NULL, NULL );
          clSafe(err, "setting signpost write to 0");
        }
      }
      forwardbackward_gpu(backward_marker);
      if(debug_rescalings){
        err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.backward_underflows_buffer, CL_TRUE, 0, 1*sizeof(int),&backward_rescalings);
        clSafe(err, "read rescalings");
        ofs<<"compute signpost marker "<<backward_marker<<", backward rescalings: "<<backward_rescalings<<endl;
      }
      if (writesignpost){
        ofs<<"for marker "<<backward_marker<<" storing sign post"<<endl;
        if (use_compressed){
          err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.unconditioned_prob_vec_buffer, CL_TRUE, 0, (compressed_states)*sizeof(packedshort_t),short_current_prob_vec);
          clSafe(err, "read current prob_vec");
          flush2disk(ofs_cache,short_current_prob_vec,compressed_states);
        }else{
          if(use_newfb){
            err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.unconditioned_prob_vec_buffer, CL_TRUE, 0, (states)*sizeof(float),current_prob_vec);
          }else{
            err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.current_prob_vec_buffer, CL_TRUE, 0, (states)*sizeof(float),current_prob_vec);
          }
          clSafe(err, "read current prob_vec");
          flush2disk(ofs_cache,current_prob_vec,states);
        }
        if(use_newfb){
          writesignpost = 0;
          err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.writesignpost_buffer, CL_TRUE, 0,  sizeof(int)*1, &writesignpost , NULL, NULL );
          clSafe(err, "setting signpost write to 0");
        }
      }
    }
    #endif
  }else{
    if(use_compressed){
      for(int cstate = 0;cstate<compressed_states;++cstate){
        short_current_prob_vec[cstate] = short_prior_vec[cstate];
      }
    }else{
      for(int state = 0;state<states;++state){
        current_prob_vec[state] = prior_vec[state];
      }
    }
    double start = clock();
    for(int marker=markers-2;marker>=(matrix_compression_ratio-1);--marker){
      //ofs<<"matrix signpost backward marker "<<marker<<endl;
      if(use_compressed){
        forwardbackward_cpu(marker,short_current_prob_vec,short_current_prob_vec,site_thetas[marker+1],backward_rescalings);
      }else{
        forwardbackward_cpu(marker,current_prob_vec,current_prob_vec,site_thetas[marker+1],backward_rescalings);
      }
      if ((marker+1) % matrix_compression_ratio == 0){
        if(use_compressed){
          flush2disk(ofs_cache,short_current_prob_vec,compressed_states);
        }else{
          flush2disk(ofs_cache,current_prob_vec,states);
        }
        cerr<<"for marker "<<marker<<" flushed to disk at row "<<signpost<<endl;
        ++signpost;
      }
    }
      // convert short to working copy of prev prob
    cerr<<"Backward rescalings in matrix signpost: "<<backward_rescalings<<endl;
    cerr<<"Backward time: "<<(clock()-start)/CLOCKS_PER_SEC<<endl;
  }
  ofs_cache.close();
  ofs<<"Closed cache output file\n";
}

void HMM_impute::init_gpu(){
  cerr<<"BLOCKWIDTH is "<<BLOCK_WIDTH<<" and SMALL_BLOCK_WIDTH is "<<SMALL_BLOCK_WIDTH<<endl;
  state_smallchunks = states/SMALL_BLOCK_WIDTH+(states%SMALL_BLOCK_WIDTH!=0);
  state_bigchunks = states/BLOCK_WIDTH+(states%BLOCK_WIDTH!=0);
  haplotype_bigchunks = total_template_haplotypes/BLOCK_WIDTH+(total_template_haplotypes%BLOCK_WIDTH!=0);
  haplotype_smallchunks = total_template_haplotypes/SMALL_BLOCK_WIDTH+(total_template_haplotypes%SMALL_BLOCK_WIDTH!=0);
  opencl_info.state_smallgroupwidth = (states/SMALL_BLOCK_WIDTH+(states%SMALL_BLOCK_WIDTH!=0))*SMALL_BLOCK_WIDTH;
  opencl_info.state_groupwidth = (states/BLOCK_WIDTH+(states%BLOCK_WIDTH!=0))*BLOCK_WIDTH;
  opencl_info.haplotypes_groupwidth = total_template_haplotypes*SMALL_BLOCK_WIDTH;
  ofs<<"Matrix signposts: "<<matrix_signposts<<" state chunks "<<state_bigchunks<<" haplotype bigchunks: "<<haplotype_bigchunks<<" state groupwidth: "<<opencl_info.state_groupwidth<<endl;
  //matrix_signposts = markers/matrix_compression_ratio+(markers%matrix_compression_ratio!=0);
  opencl_info.haplotypes_groupwidth = total_template_haplotypes*SMALL_BLOCK_WIDTH;
  ofs<<"Matrix signposts: "<<matrix_signposts<<" state chunks "<<state_bigchunks<<" haplotype bigchunks: "<<haplotype_bigchunks<<endl;
  //matrix_signposts = markers/matrix_compression_ratio+(markers%matrix_compression_ratio!=0);
  cl_int err;
  vector<cl::Platform> platforms;
  err = cl::Platform::get(&platforms);
  // Iterate over platforms
  ofs<<"Querying available OpenCL platforms: " << platforms.size() << endl;
  if (settings->platform_id>=platforms.size()) throw "Invalid platform ID";
  for(uint i=0;i<platforms.size();++i){
    ofs<<"Platform ID "<<i<<" has name "<<platforms[i].getInfo<CL_PLATFORM_NAME>().c_str()<<endl;
    if (settings->platform_id==i) ofs<<"*** SELECTED PLATFORM ABOVE ***\n";
    vector<cl::Device> devices;
    platforms[i].getDevices(CL_DEVICE_TYPE_ALL, &devices);
    ofs<<"with total devices: " << devices.size() << endl;
    //cerr<<mpi_rank<<" Querying for device "<<settings->device_id<<endl;
    if (settings->platform_id==i && settings->device_id>=devices.size()) throw "Invalid device ID";
    for (uint j = 0; j <devices.size(); ++j) {
      ofs<<"Device "<<j<<" is type ";
      string deviceName = devices[j].getInfo<CL_DEVICE_NAME>();
      cl_device_type dtype = devices[j].getInfo<CL_DEVICE_TYPE>();
      switch(dtype){
        case CL_DEVICE_TYPE_GPU:
          ofs<<"GPU\n";
          break;
        case CL_DEVICE_TYPE_CPU:
          ofs<<"CPU\n";
          break;
      }
      if (settings->platform_id==i && settings->device_id==j) ofs<<"*** SELECTED DEVICE ABOVE ***\n";
    }
  }
  ofs.flush();
   // create a context from specified platform
  cl_context_properties cps[3] = {CL_CONTEXT_PLATFORM,(cl_context_properties)(platforms[settings->platform_id])(),0};
  cl::Context context(CL_DEVICE_TYPE_GPU,cps);
  // extract the devices from this platform
  vector<cl::Device> devices = context.getInfo<CL_CONTEXT_DEVICES>();
  // make a command queue from the first device
  opencl_info.command_queue = cl::CommandQueue(context,devices[settings->device_id],0,&err);
  // Create a program from the kernel source
  ifstream ifs(settings->kernel_path.data());
  clSafe(ifs.is_open()?CL_SUCCESS:-1,"kernel_path not found");
  string source_str(istreambuf_iterator<char>(ifs),(istreambuf_iterator<char>()));
  cl::Program::Sources source(1,make_pair(source_str.c_str(),source_str.length()+1));
  cl::Program program(context,source);
  err = program.build(devices);
  if (err!=CL_SUCCESS) {
      cerr<<"Build failed:\n";
      string buffer;
      program.getBuildInfo(devices[0],CL_PROGRAM_BUILD_LOG,&buffer);
      cerr<<buffer<<endl;
  }
 // Create the OpenCL kernels
  if(use_compressed){
    opencl_info.kernel_init_forwardbackward_short = cl::Kernel(program, "init_forwardbackward_short", &err);
    clSafe(err,"create kernel init_forwardbackward_short");
    opencl_info.kernel_forward_backward_short = cl::Kernel(program, "forward_backward_short", &err);
    clSafe(err,"create kernel forward_backward");
    opencl_info.kernel_loglikelihood_short = cl::Kernel(program, "loglikelihood_short", &err);
    clSafe(err,"create kernel loglikelihood");
    opencl_info.kernel_collect_xover_marginals_block_short = cl::Kernel(program, "collect_xover_marginals_block_short", &err);
    clSafe(err,"create kernel collect_xover_marginals_block");
  }else{
    opencl_info.kernel_init_forwardbackward_float = cl::Kernel(program, "init_forwardbackward_float", &err);
    clSafe(err,"create kernel init_forwardbackward_float");
    opencl_info.kernel_forward_backward_float = cl::Kernel(program, "forward_backward_float", &err);
    clSafe(err,"create kernel forward_backward");
    opencl_info.kernel_forward_backward_emission_float = cl::Kernel(program, "forward_backward_emission_float", &err);
    clSafe(err,"create kernel forward_backward emission");
    opencl_info.kernel_loglikelihood_float = cl::Kernel(program, "loglikelihood_float", &err);
    clSafe(err,"create kernel loglikelihood");
    opencl_info.kernel_collect_xover_marginals_block_float = cl::Kernel(program, "collect_xover_marginals_block_float", &err);
    clSafe(err,"create kernel collect_xover_marginals_block");
  }
  opencl_info.kernel_reduce_fb_marginals = cl::Kernel(program, "reduce_fb_marginals", &err);
  clSafe(err,"create kernel reduce_fb_marginals");
  opencl_info.kernel_collect_xover_marginals2_block = cl::Kernel(program, "collect_xover_marginals2_block", &err);
  clSafe(err,"create kernel collect_xover_marginals_block");


// Create the memory buffers

  if (use_compressed){
    opencl_info.current_prob_vec_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (compressed_states)*sizeof(packedshort_t), NULL, &err);
    clSafe(err,"creating buffer");
    opencl_info.prior_vec_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (compressed_states)*sizeof(packedshort_t), NULL, &err);
    clSafe(err,"creating prior vec buffer");
    err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.prior_vec_buffer, CL_TRUE, 0,  sizeof(packedshort_t)*(compressed_states), short_prior_vec , NULL, NULL );
    clSafe(err, "write buffer for prior vec");
    opencl_info.backward_block_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (matrix_compression_ratio*compressed_states)*sizeof(packedshort_t), NULL, &err);
    clSafe(err,"creating backward block buffer");
  }else{
    opencl_info.current_prob_vec_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (states)*sizeof(float), NULL, &err);
    clSafe(err,"creating buffer");
    opencl_info.writesignpost_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, 1*sizeof(int), NULL, &err);
    clSafe(err,"creating buffer");
    if (use_newfb){
      opencl_info.unconditioned_prob_vec_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (states)*sizeof(float), NULL, &err);
      clSafe(err,"creating buffer");
    }
    opencl_info.prior_vec_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (states)*sizeof(float), NULL, &err);
    clSafe(err,"creating prior vec buffer");
    err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.prior_vec_buffer, CL_TRUE, 0,  sizeof(float)*(states), prior_vec , NULL, NULL );
    clSafe(err, "write buffer for prior vec");
    opencl_info.backward_block_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (matrix_compression_ratio*states)*sizeof(float), NULL, &err);
    clSafe(err,"creating backward block buffer");
  }
  opencl_info.current_subject_buffer = cl::Buffer(context, CL_MEM_READ_ONLY, (1)*sizeof(int), NULL, &err);
  clSafe(err,"creating subject buffer");
  opencl_info.marker_offset_buffer = cl::Buffer(context, CL_MEM_READ_ONLY, (1)*sizeof(int), NULL, &err);
  clSafe(err,"creating marker offset buffer");
  opencl_info.short_template_haplotypes_buffer = cl::Buffer(context, CL_MEM_READ_ONLY, (markers*total_template_haplotype_pairs)*sizeof(packedshort_t), NULL, &err);
  clSafe(err,"creating template haplo buffer");
  opencl_info.call_probs_buffer = cl::Buffer(context, CL_MEM_READ_ONLY, (all_samples*3*markers)*sizeof(float), NULL, &err);
  clSafe(err,"creating call probs buffer");
  err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.call_probs_buffer, CL_TRUE, 0,  sizeof(float)*(all_samples*3*markers), call_probs , NULL, NULL );
  clSafe(err, "write buffer for call probs");
  opencl_info.penetrance_matrix_buffer = cl::Buffer(context, CL_MEM_READ_ONLY, (markers*9)*sizeof(float), NULL, &err);
  clSafe(err,"creating pen matrix buffer");
  if (use_newfb){
    opencl_info.haplotype_marginals_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (2*total_template_haplotypes)*sizeof(float), NULL, &err);
  }else{
    opencl_info.haplotype_marginals_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (total_template_haplotypes)*sizeof(float), NULL, &err);
  }
  clSafe(err,"creating backward marginals buffer");

  opencl_info.full_marginal_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (1)*sizeof(float), NULL, &err);
  clSafe(err,"creating both marginals buffer");
  opencl_info.xover_haplotype_marginals_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (4*total_template_haplotypes)*sizeof(float), NULL, &err);
  clSafe(err,"creating xover marginals block buffer");
  opencl_info.xover_full_marginal_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (4)*sizeof(float), NULL, &err);
  clSafe(err,"creating xover marginals block buffer");
  opencl_info.marker_index_buffer = cl::Buffer(context, CL_MEM_READ_ONLY, (1)*sizeof(int), NULL, &err);
  clSafe(err,"creating backward marker buffer");
  opencl_info.emission_vec_buffer = cl::Buffer(context, CL_MEM_READ_ONLY, (states)*sizeof(float), NULL, &err);
  clSafe(err,"creating emission vec buffer");
  opencl_info.site_thetas_buffer = cl::Buffer(context, CL_MEM_READ_ONLY, (markers)*sizeof(float), NULL, &err);
  clSafe(err,"creating site thetas buffer");
  opencl_info.backward_underflows_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (1)*sizeof(int), NULL, &err);
  clSafe(err,"creating backward underflows buffer");
  opencl_info.signpost_index_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (1)*sizeof(int), NULL, &err);
  clSafe(err,"creating signpost index buffer");
  
  opencl_info.loglikelihood_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (1)*sizeof(float), NULL, &err);
  clSafe(err,"creating likelihood buffer");
  opencl_info.hap_pair_buffer = cl::Buffer(context, CL_MEM_READ_ONLY, 2*sizeof(int), NULL, &err);
  clSafe(err,"creating buffer hap pair");
  opencl_info.marginal_type_buffer = cl::Buffer(context, CL_MEM_READ_ONLY, sizeof(int), NULL, &err);
  clSafe(err,"creating buffer marginal_type");

  int arg;
  int kernelWorkGroupSize;

  if(use_compressed){
    // BEGIN SHORTKERNELS
    arg = 0;
    err = opencl_info.kernel_init_forwardbackward_short.setArg(arg++, total_template_haplotypes);
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_short.setArg(arg++, total_template_haplotype_pairs);
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_short.setArg(arg++, markers);
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_short.setArg(arg++, compressed_states);
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_short.setArg(arg++, opencl_info.current_subject_buffer);
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_short.setArg(arg++, opencl_info.marker_index_buffer);
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_short.setArg(arg++, opencl_info.marker_offset_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_init_forwardbackward_short.setArg(arg++, opencl_info.call_probs_buffer);
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_short.setArg(arg++, opencl_info.penetrance_matrix_buffer);
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_short.setArg(arg++, opencl_info.emission_vec_buffer);
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_short.setArg(arg++, opencl_info.current_prob_vec_buffer);
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_short.setArg(arg++, opencl_info.backward_block_buffer);
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_short.setArg(arg++, opencl_info.short_template_haplotypes_buffer);
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_short.setArg(arg++, opencl_info.haplotype_marginals_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_init_forwardbackward_short.setArg(arg++, cl::__local(3*sizeof(float)));
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_short.setArg(arg++, cl::__local(9*sizeof(float)));
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_short.setArg(arg++, cl::__local(BLOCK_WIDTH*sizeof(float)));
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_short.setArg(arg++, cl::__local(BLOCK_WIDTH/2*sizeof(packedshort_t)));
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    //kernelWorkGroupSize = opencl_info.kernel_init_forwardbackward_short.getWorkGroupInfo<CL_KERNEL_WORK_GROUP_SIZE>(devices[settings->device_id], &err);
    //clSafe(err,"get workgroup size kernel compute emissionvec");
    //cerr<<"Kernel compute emission vec work group size is "<<kernelWorkGroupSize<<endl;
    arg = 0;
    err = opencl_info.kernel_forward_backward_short.setArg(arg++, total_template_haplotypes);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_forward_backward_short.setArg(arg++, total_template_haplotype_pairs);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_forward_backward_short.setArg(arg++, compressed_states);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_forward_backward_short.setArg(arg++, opencl_info.marker_index_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_forward_backward_short.setArg(arg++, opencl_info.marker_offset_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_forward_backward_short.setArg(arg++, opencl_info.site_thetas_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_forward_backward_short.setArg(arg++, opencl_info.backward_underflows_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_forward_backward_short.setArg(arg++, opencl_info.emission_vec_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_forward_backward_short.setArg(arg++, opencl_info.haplotype_marginals_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_forward_backward_short.setArg(arg++, opencl_info.full_marginal_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_forward_backward_short.setArg(arg++, opencl_info.current_prob_vec_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_forward_backward_short.setArg(arg++, opencl_info.backward_block_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_forward_backward_short.setArg(arg++, cl::__local(BLOCK_WIDTH/2*sizeof(packedshort_t)));
    clSafe(err,"clSetKernelArg");
    kernelWorkGroupSize = opencl_info.kernel_forward_backward_short.getWorkGroupInfo<CL_KERNEL_WORK_GROUP_SIZE>(devices[settings->device_id], &err);
    clSafe(err,"get workgroup size kernel FB");
    //cerr<<"Kernel forwardbackward work group size is "<<kernelWorkGroupSize<<endl;
    arg = 0;
    err = opencl_info.kernel_loglikelihood_short.setArg(arg++, states);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_loglikelihood_short.setArg(arg++, opencl_info.backward_underflows_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_loglikelihood_short.setArg(arg++, opencl_info.loglikelihood_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_loglikelihood_short.setArg(arg++, opencl_info.prior_vec_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_loglikelihood_short.setArg(arg++, opencl_info.backward_block_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_loglikelihood_short.setArg(arg++, cl::__local(sizeof(packedshort_t)*BLOCK_WIDTH/2));
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_loglikelihood_short.setArg(arg++, cl::__local(sizeof(packedshort_t)*BLOCK_WIDTH/2));
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_loglikelihood_short.setArg(arg++, cl::__local(sizeof(float)*BLOCK_WIDTH));
    clSafe(err,"clSetKernelArg");
    arg = 0;
    err = opencl_info.kernel_collect_xover_marginals_block_short.setArg(arg++, compressed_states);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_collect_xover_marginals_block_short.setArg(arg++, total_template_haplotypes);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_collect_xover_marginals_block_short.setArg(arg++, total_template_haplotype_pairs);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_collect_xover_marginals_block_short.setArg(arg++, opencl_info.marker_index_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_collect_xover_marginals_block_short.setArg(arg++, opencl_info.hap_pair_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_collect_xover_marginals_block_short.setArg(arg++, opencl_info.backward_block_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_collect_xover_marginals_block_short.setArg(arg++, opencl_info.xover_haplotype_marginals_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_collect_xover_marginals_block_short.setArg(arg++, cl::__local(sizeof(packedshort_t)*BLOCK_WIDTH/2));
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_collect_xover_marginals_block_short.setArg(arg++, cl::__local(sizeof(float)*BLOCK_WIDTH));
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_collect_xover_marginals_block_short.setArg(arg++, cl::__local(sizeof(float)*BLOCK_WIDTH));
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_collect_xover_marginals_block_short.setArg(arg++, cl::__local(sizeof(float)*BLOCK_WIDTH));
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_collect_xover_marginals_block_short.setArg(arg++, cl::__local(sizeof(float)*BLOCK_WIDTH));
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_collect_xover_marginals_block_short.setArg(arg++, cl::__local(sizeof(float)*BLOCK_WIDTH));
    clSafe(err,"clSetKernelArg");
    //kernelWorkGroupSize = opencl_info.kernel_collect_xover_marginals_block.getWorkGroupInfo<CL_KERNEL_WORK_GROUP_SIZE>(devices[settings->device_id], &err);
    //clSafe(err,"get workgroup size kernel collect xover marginals");
    //ofs<<"Kernel collect xover marginals work group size is "<<kernelWorkGroupSize<<endl;
  }else{
    // BEGIN FLOATKERNELS
    arg = 0;
    err = opencl_info.kernel_init_forwardbackward_float.setArg(arg++, total_template_haplotypes);
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_float.setArg(arg++, markers);
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_float.setArg(arg++, states);
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_float.setArg(arg++, opencl_info.current_subject_buffer);
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_float.setArg(arg++, opencl_info.marker_index_buffer);
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_float.setArg(arg++, opencl_info.marker_offset_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_init_forwardbackward_float.setArg(arg++, opencl_info.call_probs_buffer);
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_float.setArg(arg++, opencl_info.penetrance_matrix_buffer);
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_float.setArg(arg++, opencl_info.emission_vec_buffer);
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_float.setArg(arg++, opencl_info.current_prob_vec_buffer);
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_float.setArg(arg++, opencl_info.backward_block_buffer);
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_float.setArg(arg++, opencl_info.short_template_haplotypes_buffer);
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_float.setArg(arg++, opencl_info.haplotype_marginals_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_init_forwardbackward_float.setArg(arg++, cl::__local(3*sizeof(float)));
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_float.setArg(arg++, cl::__local(9*sizeof(float)));
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_float.setArg(arg++, cl::__local(BLOCK_WIDTH*sizeof(float)));
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    err = opencl_info.kernel_init_forwardbackward_float.setArg(arg++, cl::__local(BLOCK_WIDTH/2*sizeof(packedshort_t)));
    clSafe(err,"clSetKernelArgComputeEmissionVector2");
    kernelWorkGroupSize = opencl_info.kernel_init_forwardbackward_float.getWorkGroupInfo<CL_KERNEL_WORK_GROUP_SIZE>(devices[settings->device_id], &err);
    clSafe(err,"get workgroup size kernel compute emissionvec");
    ofs<<"Kernel compute emission vec work group size is "<<kernelWorkGroupSize<<endl;

    arg = 0;
    err = opencl_info.kernel_forward_backward_emission_float.setArg(arg++, total_template_haplotypes);
    clSafe(err,"clSetKernelArg_forwardbackwardemission");
    err = opencl_info.kernel_forward_backward_emission_float.setArg(arg++, markers);
    clSafe(err,"clSetKernelArg_forwardbackwardemission");
    err = opencl_info.kernel_forward_backward_emission_float.setArg(arg++, states);
    clSafe(err,"clSetKernelArg_forwardbackwardemission");
    err = opencl_info.kernel_forward_backward_emission_float.setArg(arg++, opencl_info.current_subject_buffer);
    clSafe(err,"clSetKernelArg_forwardbackwardemission");
    err = opencl_info.kernel_forward_backward_emission_float.setArg(arg++, opencl_info.marker_index_buffer);
    clSafe(err,"clSetKernelArg_forwardbackwardemission");
    err = opencl_info.kernel_forward_backward_emission_float.setArg(arg++, opencl_info.marker_offset_buffer);
    clSafe(err,"clSetKernelArg_forwardbackwardemission");
    err = opencl_info.kernel_forward_backward_emission_float.setArg(arg++, opencl_info.writesignpost_buffer);
    clSafe(err,"clSetKernelArg_forwardbackwardemission");
    err = opencl_info.kernel_forward_backward_emission_float.setArg(arg++, opencl_info.site_thetas_buffer);
    clSafe(err,"clSetKernelArg_forwardbackwardemission");
    err = opencl_info.kernel_forward_backward_emission_float.setArg(arg++, opencl_info.backward_underflows_buffer);
    clSafe(err,"clSetKernelArg_forwardbackwardemission");
    err = opencl_info.kernel_forward_backward_emission_float.setArg(arg++, opencl_info.call_probs_buffer);
    clSafe(err,"clSetKernelArg_forwardbackwardemission");
    err = opencl_info.kernel_forward_backward_emission_float.setArg(arg++, opencl_info.penetrance_matrix_buffer);
    clSafe(err,"clSetKernelArg_forwardbackwardemission");
    err = opencl_info.kernel_forward_backward_emission_float.setArg(arg++, opencl_info.short_template_haplotypes_buffer);
    clSafe(err,"clSetKernelArg_forwardbackwardemission");
    err = opencl_info.kernel_forward_backward_emission_float.setArg(arg++, opencl_info.emission_vec_buffer);
    clSafe(err,"clSetKernelArg_forwardbackwardemission");
    err = opencl_info.kernel_forward_backward_emission_float.setArg(arg++, opencl_info.haplotype_marginals_buffer);
    clSafe(err,"clSetKernelArg_forwardbackwardemission");
    err = opencl_info.kernel_forward_backward_emission_float.setArg(arg++, opencl_info.full_marginal_buffer);
    clSafe(err,"clSetKernelArg_forwardbackwardemission");
    err = opencl_info.kernel_forward_backward_emission_float.setArg(arg++, opencl_info.current_prob_vec_buffer);
    clSafe(err,"clSetKernelArg_forwardbackwardemission");
    err = opencl_info.kernel_forward_backward_emission_float.setArg(arg++, opencl_info.backward_block_buffer);
    clSafe(err,"clSetKernelArg_forwardbackwardemission");
    err = opencl_info.kernel_forward_backward_emission_float.setArg(arg++, opencl_info.unconditioned_prob_vec_buffer);
    clSafe(err,"clSetKernelArg_forwardbackwardemission");
    err = opencl_info.kernel_forward_backward_emission_float.setArg(arg++, cl::__local(3*sizeof(float)));
    clSafe(err,"clSetKernelArg_forwardbackwardemission");
    err = opencl_info.kernel_forward_backward_emission_float.setArg(arg++, cl::__local(9*sizeof(float)));
    clSafe(err,"clSetKernelArg_forwardbackwardemission");
    err = opencl_info.kernel_forward_backward_emission_float.setArg(arg++, cl::__local(BLOCK_WIDTH*sizeof(float)));
    clSafe(err,"clSetKernelArg_forwardbackwardemission");
    err = opencl_info.kernel_forward_backward_emission_float.setArg(arg++, cl::__local(BLOCK_WIDTH*sizeof(float)));
    clSafe(err,"clSetKernelArg_forwardbackwardemission");
    err = opencl_info.kernel_forward_backward_emission_float.setArg(arg++, cl::__local(BLOCK_WIDTH/2*sizeof(packedshort_t)));
    clSafe(err,"clSetKernelArg_forwardbackwardemission");
    kernelWorkGroupSize = opencl_info.kernel_forward_backward_emission_float.getWorkGroupInfo<CL_KERNEL_WORK_GROUP_SIZE>(devices[settings->device_id], &err);
    clSafe(err,"get workgroup size kernel FB");
    ofs<<"Kernel forwardbackward emission work group size is "<<kernelWorkGroupSize<<endl;
    arg = 0;
    err = opencl_info.kernel_forward_backward_float.setArg(arg++, total_template_haplotypes);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_forward_backward_float.setArg(arg++, states);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_forward_backward_float.setArg(arg++, opencl_info.marker_index_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_forward_backward_float.setArg(arg++, opencl_info.marker_offset_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_forward_backward_float.setArg(arg++, opencl_info.site_thetas_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_forward_backward_float.setArg(arg++, opencl_info.backward_underflows_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_forward_backward_float.setArg(arg++, opencl_info.emission_vec_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_forward_backward_float.setArg(arg++, opencl_info.haplotype_marginals_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_forward_backward_float.setArg(arg++, opencl_info.full_marginal_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_forward_backward_float.setArg(arg++, opencl_info.current_prob_vec_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_forward_backward_float.setArg(arg++, opencl_info.backward_block_buffer);
    clSafe(err,"clSetKernelArg");
    kernelWorkGroupSize = opencl_info.kernel_forward_backward_float.getWorkGroupInfo<CL_KERNEL_WORK_GROUP_SIZE>(devices[settings->device_id], &err);
    clSafe(err,"get workgroup size kernel FB");
    ofs<<"Kernel forwardbackward work group size is "<<kernelWorkGroupSize<<endl;


    arg = 0;
    err = opencl_info.kernel_loglikelihood_float.setArg(arg++, states);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_loglikelihood_float.setArg(arg++, opencl_info.backward_underflows_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_loglikelihood_float.setArg(arg++, opencl_info.loglikelihood_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_loglikelihood_float.setArg(arg++, opencl_info.prior_vec_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_loglikelihood_float.setArg(arg++, opencl_info.backward_block_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_loglikelihood_float.setArg(arg++, cl::__local(sizeof(float)*BLOCK_WIDTH));
    clSafe(err,"clSetKernelArg");

    arg = 0;
    err = opencl_info.kernel_collect_xover_marginals_block_float.setArg(arg++, states);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_collect_xover_marginals_block_float.setArg(arg++, total_template_haplotypes);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_collect_xover_marginals_block_float.setArg(arg++, opencl_info.marker_index_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_collect_xover_marginals_block_float.setArg(arg++, opencl_info.hap_pair_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_collect_xover_marginals_block_float.setArg(arg++, opencl_info.backward_block_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_collect_xover_marginals_block_float.setArg(arg++, opencl_info.xover_haplotype_marginals_buffer);
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_collect_xover_marginals_block_float.setArg(arg++, cl::__local(sizeof(float)*BLOCK_WIDTH));
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_collect_xover_marginals_block_float.setArg(arg++, cl::__local(sizeof(float)*BLOCK_WIDTH));
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_collect_xover_marginals_block_float.setArg(arg++, cl::__local(sizeof(float)*BLOCK_WIDTH));
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_collect_xover_marginals_block_float.setArg(arg++, cl::__local(sizeof(float)*BLOCK_WIDTH));
    clSafe(err,"clSetKernelArg");
    err = opencl_info.kernel_collect_xover_marginals_block_float.setArg(arg++, cl::__local(sizeof(float)*BLOCK_WIDTH));
    clSafe(err,"clSetKernelArg");
    //kernelWorkGroupSize = opencl_info.kernel_collect_xover_marginals_block.getWorkGroupInfo<CL_KERNEL_WORK_GROUP_SIZE>(devices[settings->device_id], &err);
    //clSafe(err,"get workgroup size kernel collect xover marginals");
    //ofs<<"Kernel collect xover marginals work group size is "<<kernelWorkGroupSize<<endl;
  }

  arg = 0;
  err = opencl_info.kernel_reduce_fb_marginals.setArg(arg++, total_template_haplotypes);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_reduce_fb_marginals.setArg(arg++, opencl_info.marker_index_buffer);
    clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_reduce_fb_marginals.setArg(arg++, opencl_info.haplotype_marginals_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_reduce_fb_marginals.setArg(arg++, opencl_info.full_marginal_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_reduce_fb_marginals.setArg(arg++, cl::__local(sizeof(float) * BLOCK_WIDTH ));
  clSafe(err,"clSetKernelArg");
  kernelWorkGroupSize = opencl_info.kernel_reduce_fb_marginals.getWorkGroupInfo<CL_KERNEL_WORK_GROUP_SIZE>(devices[settings->device_id], &err);
  clSafe(err,"get workgroup size kernel reduce fb marginals second");
  //cerr<<"Kernel reduce 2nd marginals work group size is "<<kernelWorkGroupSize<<endl;

  arg = 0;
  err = opencl_info.kernel_collect_xover_marginals2_block.setArg(arg++, total_template_haplotypes);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_collect_xover_marginals2_block.setArg(arg++, opencl_info.xover_haplotype_marginals_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_collect_xover_marginals2_block.setArg(arg++, opencl_info.xover_full_marginal_buffer);
  clSafe(err,"clSetKernelArg");
  //err = opencl_info.kernel_collect_xover_marginals2_block.setArg(arg++, cl::__local(sizeof(float)*SMALL_BLOCK_WIDTH));
  err = opencl_info.kernel_collect_xover_marginals2_block.setArg(arg++, cl::__local(sizeof(float)*BLOCK_WIDTH));
  clSafe(err,"clSetKernelArg");
  //kernelWorkGroupSize = opencl_info.kernel_collect_xover_marginals2_block.getWorkGroupInfo<CL_KERNEL_WORK_GROUP_SIZE>(devices[settings->device_id], &err);
  //clSafe(err,"get workgroup size kernel collect xover marginals");
  //ofs<<"Kernel collect xover marginals work group size is "<<kernelWorkGroupSize<<endl;
  //exit(0);
}

int main_hmm_impute(int argc, char * argv[]){
  return 0;
}
